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

Question

我有一个一维的数据数组，想要提取它的空间变化。通常的做法是对数据进行移动线性回归，然后保存它的斜率...

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

我之前用过一些来自比较老的scikits.timeseries模块的移动窗口函数，效果还不错。这些函数是用C语言写的，但我还没能在窗口大小变化的情况下使用它们（不确定你是否需要这个功能）。

http://pytseries.sourceforge.net/lib.moving_funcs.html

如果你需要下载，可以去这里（如果你用的是Python 2.7或更高版本，可能需要自己编译这个扩展——我在2.7上做过，效果很好）：

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

如果你能稍微整理一下你的示例代码，我们可能能更好地帮助你。我建议你在第7和第8行定义一些参数/对象（比如你定义的'small'）时，把它们当作变量，这样第8行就不会有那么多难以理解的括号了。

Answer 2

好的，我找到了一种看起来像是解决办法的方法……虽然不算是直接回答问题，但这是一个进行多点移动微分的方式。我测试过这个方法，结果看起来和移动回归非常相似……我使用了一个一维的索贝尔滤波器（就是一个从-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

线性回归的计算是基于各种值的总和。所以你可以写一个更高效的程序，在窗口移动时更新这个总和（也就是加上一个新点，同时减去一个早先的点）。

这样做比每次窗口移动时都重新计算要高效得多，但可能会出现四舍五入的误差。因此，你需要偶尔重新开始计算。

如果数据点间隔相等，你可能可以通过提前计算所有的x相关值来做得更好，不过我对你的例子不太了解，所以不确定这是否适用。

所以我想我就假设它是适用的。

斜率的计算公式是 (NΣXY - (ΣX)(ΣY)) / (NΣX² - (ΣX)²)，这里的“²”表示平方 - 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_N“进入”窗口。