一次性对 scipy 的“curve_fit”进行多次迭代

2023-12-07

考虑以下 MWE

import numpy as np
from scipy.optimize import curve_fit
X=np.arange(1,10,1)
Y=abs(X+np.random.randn(15,9))

def linear(x, a, b):
    return (x/b)**a

coeffs=[]
for ix in range(Y.shape[0]):
    print(ix)
    c0, pcov = curve_fit(linear, X, Y[ix])
    coeffs.append(c0)


XX=np.tile(X, Y.shape[0])
c0, pcov = curve_fit(linear, XX, Y.flatten())

我遇到一个问题,我基本上必须这样做,但不是 15 次迭代,而是数千次迭代,而且速度非常慢。

有没有什么办法可以同时完成所有这些迭代curve_fit?我知道函数的结果应该是一维数组,所以只需像这样传递参数

c0, pcov = curve_fit(nlinear, X, Y)

是行不通的。我也认为答案必须是扁平化Y,所以我可以获得平坦的结果,但我就是无法让任何东西发挥作用。

EDIT

我知道如果我做类似的事情

XX=np.tile(X, Y.shape[0])
c0, pcov = curve_fit(nlinear, XX, Y.flatten())

然后我得到系数的“平均值”,但这不是我想要的。

EDIT 2

作为记录,我使用 Jacques Kvam 的设置解决了问题,但使用 Numpy 实现了(由于限制)

lX = np.log(X)
lY = np.log(Y)
A = np.vstack([lX, np.ones(len(lX))]).T
m, c=np.linalg.lstsq(A, lY.T)[0]

进而m is a并得到b:

b=np.exp(-c/m)

最小二乘不会给出相同的结果,因为在这种情况下噪声是通过对数转换的。如果噪声为零,则两种方法都会给出相同的结果。

import numpy as np
from numpy import random as rng
from scipy.optimize import curve_fit
rng.seed(0)
X=np.arange(1,7)
Y = np.zeros((4, 6))
for i in range(4):
    b = a = i + 1
    Y[i] = (X/b)**a + 0.01 * randn(6)

def linear(x, a, b):
    return (x/b)**a

coeffs=[]
for ix in range(Y.shape[0]):
    print(ix)
    c0, pcov = curve_fit(linear, X, Y[ix])
    coeffs.append(c0)

coefs is

[array([ 0.99309127,  0.98742861]),
 array([ 2.00197613,  2.00082722]),
 array([ 2.99130237,  2.99390585]),
 array([ 3.99644048,  3.9992937 ])]

我将使用 scikit-learn 的线性回归实现,因为我相信它可以很好地扩展。

from sklearn.linear_model import LinearRegression
lr = LinearRegression()

获取日志X and Y

lX = np.log(X)[None, :]
lY = np.log(Y)

现在拟合并检查系数是否与之前相同。

lr.fit(lX.T, lY.T)
lr.coef_

这给出了相似的指数。

array([ 0.98613517,  1.98643974,  2.96602892,  4.01718514])

现在检查除数。

np.exp(-lr.intercept_ / lr.coef_.ravel())

这给出了相似的系数,您可以看到这些方法在答案中略有不同。

array([ 0.99199406,  1.98234916,  2.90677142,  3.73416501])
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

一次性对 scipy 的“curve_fit”进行多次迭代 的相关文章

随机推荐