Numpy Polyfit 或 X 和 Y 多维数组的任何拟合

Posted

技术标签:

【中文标题】Numpy Polyfit 或 X 和 Y 多维数组的任何拟合【英文标题】:Numpy Polyfit or any fitting to X and Y multidimensional arrays 【发布时间】:2014-04-25 20:37:12 【问题描述】:

我有两个大的多维数组:Y 携带三百万个对象的测量值(例如 shape=(500000,3)),X 具有相同的形状,但包含 Y 测量值的位置。

首先,我希望包含一个对象的每一行都拟合一个多项式方程。我知道遍历数组很慢,但我现在正在做的是:

fit = array([polyfit(X[i],Y[i],deg) for i in xrange(obs.shape[0])])

我的问题是:是否有可能在不显式迭代它们的情况下拟合两个数组的每一行?

【问题讨论】:

您是否尝试过沿轴应用之类的方法? 3 个点只能拟合线性或二次(精确拟合)函数。在这两种情况下都可以计算出显式​​解,然后可以对其进行矢量化以同时计算所有行。 其实三点只是举例,当我真正应用该方法时,尺寸会更大。 【参考方案1】:

可以在不沿第一个轴迭代的情况下这样做。但是,您的第二个轴相当短(只有 3 个),您实际上可以拟合不超过 2 个系数。

In [67]:

import numpy as np
import scipy.optimize as so

In [68]:

def MD_ployError(p, x, y):
    '''if x has the shape of (n,m), y must be (n,m), p must be (n*p, ), where p is degree'''
    #d is no. of degree
    p_rshp=p.reshape((x.shape[0], -1))
    f=y*1.
    for i in range(p_rshp.shape[1]):
        f-=p_rshp[:,i][:,np.newaxis]*(x**i)
    return (f**2).sum()

In [69]:

X=np.random.random((100, 6))
Y=4+2*X+3*X*X
P=(np.zeros((100,3))+[1,1,1]).ravel()

In [70]:

MD_ployError(P, X, Y)

Out[70]:
11012.2067606684

In [71]:

R=so.fmin_slsqp(MD_ployError, P, args=(X, Y))
Iteration limit exceeded    (Exit mode 9) #you can increase iteration limit, but the result is already good enough.
            Current function value: 0.00243784856039
            Iterations: 101
            Function evaluations: 30590
            Gradient evaluations: 101

In [72]:

R.reshape((100, -1))

Out[72]:
array([[ 3.94488512,  2.25402422,  2.74773571],
       [ 4.00474864,  1.97966551,  3.02010015],
       [ 3.99919559,  2.0032741 ,  2.99753804],
..............................................)

【讨论】:

谢谢朱。但是,这种方法比我的单行实现更快吗? 它依赖。如果您可以提供一个函数来计算误差函数的导数,它会更快。如果初始猜测参数是“好猜测”,它会更快。不同的优化器可能会更快。但总的来说,我认为这可能是一种矫枉过正,特别是如果你没有设置任何约束(比如 x^2 的系数必须相等)。我也遇到过同样的情况,我选择multiprocessing(但在这种情况下,我为每一行拟合了一个非常复杂的模型)。毕竟,任何一行的拟合都完全独立于其他行。【参考方案2】:

是的,如果您使用来自np.polynomial、not the old np.polyfit 的新 numpy polyfit:

X = np.arange(3)
Y = np.random.rand(10000, 3)

fit = np.array([np.polyfit(X, y, 2) for y in Y])
fits = np.polynomial.polynomial.polyfit(X, Y.T, 2)

assert np.allclose(fit.T[::-1], fits)

时间:

In [692]: timeit fit = np.array([np.polyfit(X, y, 2) for y in Y])
 1 loops, best of 3: 2.22 s per loop

In [693]:  timeit fits = np.polynomial.polynomial.polyfit(X, Y.T, 2)
100 loops, best of 3: 3.63 ms per loop

【讨论】:

我认为 OP 意味着 X 也具有 (10000,3) 的形状,对于 d 度数(系数),结果需要为 (10000,d):适合 x 的每一行由y的每一行@ 其实X数组也必须是多维的。正如我上面评论的那样,我确实想要更通用的解决方案。选择维度 3 只是为了提出问题。 维度3 不是这里的限制,它是数组的形状。也就是说,y 中的每一行都有不同的x 值。一组位置 (x) 通常有多组测量值 (y),但您拥有完全独立的测量值和位置集。我认为没有一种简单的方法可以使用内置多项式拟合一次拟合它们,您必须像 @CT 那样定义自己的误差函数。

以上是关于Numpy Polyfit 或 X 和 Y 多维数组的任何拟合的主要内容,如果未能解决你的问题,请参考以下文章

在 numpy polyfit 中使用的权重值是多少,拟合的误差是多少

Numpy:连接多维和一维数组

numpy中多维数组的绝对索引

numpy-np.ceil,np.floor,np.expand_dims方法

numpy中是否有多维版本的arange / linspace?

numpy.polyfit 和 scipy.polyfit 有啥区别? [复制]