首页 文章

线性插值 - 制作网格

提问于
浏览
1

我想在不同的模型之间进行插值 . 为方便起见,我的数据如下所示:

enter image description here

我有10种不同的模拟(我称之为 z ) . 对于每个 z ,我有一个 array x 和一个 array y (对于给定的 zlen(x)=len(y) ) . 例如:

for z=1x.shape=(1200,)y.shape=(1200,)

for z=2x.shape=(1250,)y.shape=(1250,)

for z=3x.shape=(1236,)y.shape=(1236,)

等等 ...

我想插值,以便对于给定的 zx ,我得到 y . 例如,对于 z=2.5x=10**9 ,代码输出 y . 我假设:

y = a*x + b*z + c 当然我不知道 abc .

我的问题是如何将数据存储在网格中?我很困惑,因为 x 的大小与 xy 不同 . 如何构建网格?

更新

我能够部分解决我的问题 . 我首先做的是使用 interp1dxy 之间进行插值 . 它工作得非常好 . 然后我创建了一个 xy 值的新网格 . 简而言之,该方法是:

f = interp1d(x, y, kind='linear')
new_x = np.linspace(10**7, 4*10**9, 10000)
new_y = f(new_x)

然后我插入 xyz

ff = LinearNDInterpolator( (x, z), y)

要测试方法是否有效,这是一个带有 z=3 的图 .

enter image description here

情节看起来很好,直到 x=10**8 . 实际上,该线偏离了原始模型 . 这是我进一步放大时的情节:

enter image description here

x > 10**8 时插值明显不好 . 我该如何解决?

2 回答

  • 0

    似乎在您的问题中曲线y(x)表现良好,因此您可能只是先为z的给定值插入y(x),然后在获得的y值之间进行插值 .

    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.widgets import Slider
    
    import random
    
    #####
    # Generate some data
    #####
    generate = lambda x, z: 1./(x+1.)+(z*x/75.+z/25.)
    
    def f(z):
        #create an array of values between zero and 100 of random length
        x = np.linspace(0,10., num=random.randint(42,145))
        #generate corresponding y values
        y = generate(x, z)
        return np.array([x,y])
    
    Z = [1, 2, 3, 3.6476, 4, 5.1]
    A = [f(z) for z in Z]
    #now A contains the dataset of [x,y] pairs for each z value
    
    #####
    # Interpolation
    #####
    def do_interpolation(x,z):
        #assume Z being sorted in ascending order
        #look for indizes of z values closest to given z
        ig = np.searchsorted(Z, z)
        il = ig-1
        #interpolate y(x) for those z values
        yg = np.interp(x, A[ig][0,:], A[ig][1,:])
        yl = np.interp(x, A[il][0,:], A[il][1,:])
        #linearly interpolate between yg and yl  
        return yl + (yg-yl)*float(z-Z[il])/(Z[ig] - Z[il])  
    
    # do_interpolation(x,z) will now provide the interpolated data
    print do_interpolation( np.linspace(0, 10), 2.5) 
    
    #####
    # Plotting, use Slider to change the value of z. 
    #####
    fig=plt.figure()
    fig.subplots_adjust(bottom=0.2)
    ax=fig.add_subplot(111)
    for i in range(len(Z)):
        ax.plot(A[i][0,:] , A[i][1,:], label="{z}".format(z=Z[i]) )
    
    l, = ax.plot(np.linspace(0, 10) , do_interpolation( np.linspace(0, 10), 2.5), label="{z}".format(z="interpol"), linewidth=2., color="k" )
    
    axn1 = plt.axes([0.25, 0.1, 0.65, 0.03], axisbg='#e4e4e4')
    sn1 = Slider(axn1, 'z', Z[0], Z[-1], valinit=2.5)
    def update(val):
        l.set_data(np.linspace(0, 10), do_interpolation( np.linspace(0, 10), val))
        plt.draw()
    sn1.on_changed(update)
    
    ax.legend()
    plt.show()
    
  • 0

    你正在做什么对我来说似乎有点奇怪,至少你似乎使用一组 y 值来进行插值 . 我建议不要一个接一个地执行两个插值,但考虑到你的 y(z,x) 函数是纯二维插值问题的结果 .

    正如我在评论中指出的那样,我建议使用 scipy.interpolate.LinearNDInterpolator ,这是 griddata 在引擎盖下用于双线性插值的同一个对象 . 正如我们在评论中所讨论的那样,你需要有一个插值器,你可以在之后多次查询,所以我们必须使用低级插值器对象,因为它是可调用的 .

    以下是我的意思的完整示例,包括虚拟数据和绘图:

    import numpy as np
    import scipy.interpolate as interp
    import matplotlib.pyplot as plt
    
    # create dummy data
    zlist = range(4)  # z values
    # one pair of arrays for each z value in a list:
    xlist = [np.linspace(-1,1,41),
             np.linspace(-1,1,61),
             np.linspace(-1,1,55),
             np.linspace(-1,1,51)]
    funlist = [lambda x:0.1*np.ones_like(x),
               lambda x:0.2*np.cos(np.pi*x)+0.4,
               lambda x:np.exp(-2*x**2)+0.5,
               lambda x:-0.7*np.abs(x)+1.7]
    ylist = [f(x) for f,x in zip(funlist,xlist)]
    
    # create contiguous 1d arrays for interpolation
    all_x = np.concatenate(xlist)
    all_y = np.concatenate(ylist)
    all_z = np.concatenate([np.ones_like(x)*z for x,z in zip(xlist,zlist)])
    
    # create a single linear interpolator object
    yfun = interp.LinearNDInterpolator((all_z,all_x),all_y)
    
    # generate three interpolated sets: one with z=2 to reproduce existing data,
    # two with z=1.5 and z=2.5 respectively to see what happens
    xplot = np.linspace(-1,1,30)
    z = 2
    y_repro = yfun(z,xplot)
    z = 1.5
    y_interp1 = yfun(z,xplot)
    z = 2.5
    y_interp2 = yfun(z,xplot)
    
    # plot the raw data (markers) and the two interpolators (lines)
    fig,ax = plt.subplots()
    for x,y,z,mark in zip(xlist,ylist,zlist,['s','o','v','<','^','*']):
        ax.plot(x,y,'--',marker=mark,label='z={}'.format(z))
    ax.plot(xplot,y_repro,'-',label='z=2 interp')
    ax.plot(xplot,y_interp1,'-',label='z=1.5 interp')
    ax.plot(xplot,y_interp2,'-',label='z=2.5 interp')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    # reduce plot size and put legend outside for prettiness, see also http://stackoverflow.com/a/4701285/5067311
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    
    plt.show()
    

    result

    您没有指定如何存储 (x,y) 数组对的系列,我使用了numpy ndarray 的列表 . 如您所见,我将1d数组列表压缩为一组1d数组: all_xall_yall_z . 这些可以用作分散的数据,您可以从中构造插值器对象 . 正如您在结果中看到的那样,对于 z=2 它会再现输入点,对于非整数 z ,它会在相关的 y(x) 曲线之间进行插值 .

    此方法应适用于您的数据集 . 但需要注意的是: x 轴上的对数刻度上有大量数字 . 仅这一点就可能导致数字不稳定 . 我建议您也尝试使用 log(x) 执行插值,它可能表现得更好(这只是一个模糊的猜测) .

相关问题