在Python中对一维数据进行移动线性拟合

Question

我有一个一维数据数组，想要提取空间变化。我希望将标准方法python化，即对数据执行移动线性回归并保存梯度...

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中实现移动拟合的方法?

- Scott Collis

---

不要想象，进行剖析！使用cProfile。如果这还不够细致，请尝试将部分代码拆分成自己的函数。 - Nick ODell

你试过不用打印语句吗？它们可能会使程序变慢...而且，三重嵌套的for循环会很慢。 - tylerthemiler

使用 range 而不是 xrange 可能会导致速度变慢。 - Xavier Combelle

这段代码真的能工作吗？kdp[swn, rayn, :]看起来不对，也许你的意思是kdp[swn, rayn:]？此外，变量swn和ray在实际计算中都没有被使用，这样做对吗？最后一部分使用itertools.chain看起来会产生语法错误，因为您在*[...]之后使用了变量。最后只有phidp.shape被使用，而不是phidp本身的数据——这代表什么？ - Klohkwherk

抱歉，刚刚有一些粘贴错误……并且添加了打印语句，以及做了一点小修改……大部分时间都花在 [for in] 循环上了。 - Scott Collis

哦，谢谢...我需要做更多的分析，你说得完全正确。 - Scott Collis

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_N则“进入”。

- andrew cooke

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)

这个运行速度很快！

- Scott Collis

Answer 3

我曾经在旧版本的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行使用太多难以跟踪的括号了。

- mdscruggs

---

我已经清理了代码，忽略itertools那一行...关键在于small=[polyfit(distance[a:a+2fitlen], phidp[swn, rayn, a:a+2fitlen],1)[0] for a in xrange(myshape[2]-2*fitlen)]。看起来这些函数非常接近！我在考虑使用1D Sobel滤波器... - Scott Collis