在Python中对1D数据执行移动线性拟合

Question

我有一个一维数组的数据，希望提取空间变化。我希望pythonize的标准方法是对数据执行移动线性回归并保存梯度。。。在

def nssl_kdp(phidp, distance, fitlen):
    kdp=zeros(phidp.shape, dtype=float)
    myshape=kdp.shape
    for swn in range(myshape[0]):
        print "Sweep ", swn+1
        for rayn in range(myshape[1]):
            print "ray ", rayn+1
            small=[polyfit(distance[a:a+2*fitlen], phidp[swn, rayn, a:a+2*fitlen],1)[0] for a in xrange(myshape[2]-2*fitlen)]
            kdp[swn, rayn, :]=array((list(itertools.chain(*[fitlen*[small[0]], small, fitlen*[small[-1]]]))))
    return kdp

这很有效，但速度很慢。。。我需要做17*360次。。。在

我想开销在[for in arange]行的迭代器中。。。在numpy/scipy中是否有一个移动配合的实现？在

Answer 1

线性回归的计算是以各种值之和为基础的。因此，您可以编写一个更高效的例程，在窗口移动时修改总和（加一个点，减去一个较早的点）。在

这将比每次窗口移动时重复该过程要高效得多，但可能会出现舍入错误。所以你需要偶尔重启。在

对于等间距的点，通过预先计算所有的x依赖项，您可能可以做得更好，但是我不了解您的示例的详细信息，因此不确定它是否相关。在

所以我想我会假设是的。在

斜率是（N∑XY—∑X）（∑Y））/（N∑X2—∑X）2），其中“2”是“平方”-http://easycalculation.com/statistics/learn-regression.php

对于均匀分布的数据，分母是固定的（因为您可以将x轴移动到窗口的开头而不更改渐变）。分子中的（∑X）也是固定的（出于同样的原因）。所以你只需要关心∑XY和∑Y，后者很简单，只需加减一个值。前者每一步减少∑Y（每个X权重减少1），增加（N-1）Y（假设X_0为0，X_N为N-1）。在

我想这还不清楚。我要说的是，斜率的公式不需要每一步都完全重新计算。特别是因为，在每个步骤中，可以将X值重命名为0,1，…N-1，而不更改坡度。所以公式中几乎所有的东西都是一样的。所有这些变化都是两个条件，这取决于Y，因为Y_0“退出”窗口，而Y un“移入”。在

Answer 2

好吧。。我有一个解决办法。。不是一个答案，而是一种移动的，多点微分的方法。。。我已经测试了这个结果看起来非常类似于一个移动的回归。。。我使用了1D sobel过滤器（从-1到1的斜坡与数据卷积）：

def KDP(phidp, dx, fitlen):
    kdp=np.zeros(phidp.shape, dtype=float)
    myshape=kdp.shape
    for swn in range(myshape[0]):
        #print "Sweep ", swn+1
        for rayn in range(myshape[1]):
            #print "ray ", rayn+1
            kdp[swn, rayn, :]=sobel(phidp[swn, rayn,:], window_len=fitlen)/dx
    return kdp

def sobel(x,window_len=11):
    """Sobel differential filter for calculating KDP
    output:
        differential signal (Unscaled for gate spacing
        example:
    """
    s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
    #print(len(s))
    w=2.0*np.arange(window_len)/(window_len-1.0) -1.0
    #print w
    w=w/(abs(w).sum())
    y=np.convolve(w,s,mode='valid')
    return -1.0*y[window_len/2:len(x)+window_len/2]/(window_len/3.0)

这跑得很快！在

Answer 3

我用过这些移动窗口函数scikits.timeseries模块取得了一些成功。它们是用C语言实现的，但我还没有设法在移动窗口大小不同的情况下使用它们（不确定是否需要该功能）。在

{a1}

点击这里下载（如果使用Python2.7+，您可能需要编译扩展本身我为2.7做了这个，它工作得很好）：

http://sourceforge.net/projects/pytseries/files/scikits.timeseries/0.91.3/

如果您稍微清理一下您的示例代码，我/我们可能会帮助您更多。我会考虑将第7行和第8行中的一些参数/对象定义为变量（在这里定义“small”），这样就不会在第8行末尾加上太多难以跟随的括号。在