考虑以下 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)
不会工作。另外我认为答案必须是 flattening Y
,所以我可以得到一个 flattened 结果,但我无法得到任何工作。
编辑
我知道如果我做类似的事情
XX=np.tile(X, Y.shape[0])
c0, pcov = curve_fit(nlinear, XX, Y.flatten())
然后我得到系数的“平均”值,但这不是我想要的。
编辑 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
是a
和得到b
:
b=np.exp(-c/m)
最小二乘法不会给出相同的结果,因为在这种情况下噪声是由 log 转换的。如果噪声为零,则两种方法都会给出相同的结果。
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
是
[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
和的日志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])
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句