我试图从我一直在运行的仿真代码中拟合一些数据,以找出幂律的相关性。当我绘制线性拟合时,数据拟合得不是很好。
这是我用来拟合数据的python脚本:
#!/usr/bin/env python
from scipy import optimize
import numpy
xdata=[ 0.00010851, 0.00021701, 0.00043403, 0.00086806, 0.00173611, 0.00347222]
ydata=[ 29.56241016, 29.82245508, 25.33930469, 19.97075977, 12.61276074, 7.12695312]
fitfunc = lambda p, x: p[0] + p[1] * x ** (p[2])
errfunc = lambda p, x, y: (y - fitfunc(p, x))
out,success = optimize.leastsq(errfunc, [1,-1,-0.5],args=(xdata, ydata),maxfev=3000)
print "%g + %g*x^%g"%(out[0],out[1],out[2])
我得到的输出是:
-71205.3 + 71174.5 * x ^ -9.79038e-05
虽然在情节上,拟合看起来像您从最小二乘拟合中所期望的一样好,但是输出的形式使我感到困扰。我希望常数会接近您期望的零值(大约30)。而且我期望找到比10 ^ -5大的功率相关性。
我试过重新缩放数据并使用参数来优化optimize.leastsq,但运气不佳。我要完成的任务是否可行?或者我的数据是否不允许这样做?计算成本很高,因此获取更多数据点并非易事。
谢谢!
最佳答案
它有助于重新调整xdata
的大小,因此数字并非都那么小。
您可以使用新的变量xprime = 1000*x
。
然后将xprime
与y
匹配。
最小二乘会找到q
拟合参数
y = q[0] + q[1] * (xprime ** q[2])
= q[0] + q[1] * ((1000*x) ** q[2])
所以让
p[0] = q[0]
p[1] = q[1] * (1000**q[2])
p[2] = q[2]
然后
y = p[0] + p[1] * (x ** p[2])
它还有助于将初始猜测更改为更接近您期望的结果,例如
[max(ydata), -1, -0.5]
。from scipy import optimize
import numpy as np
def fitfunc(p, x):
return p[0] + p[1] * (x ** p[2])
def errfunc(p, x, y):
return y - fitfunc(p, x)
xdata=np.array([ 0.00010851, 0.00021701, 0.00043403, 0.00086806,
0.00173611, 0.00347222])
ydata=np.array([ 29.56241016, 29.82245508, 25.33930469, 19.97075977,
12.61276074, 7.12695312])
N = 5000
xprime = xdata * N
qout,success = optimize.leastsq(errfunc, [max(ydata),-1,-0.5],
args=(xprime, ydata),maxfev=3000)
out = qout[:]
out[0] = qout[0]
out[1] = qout[1] * (N**qout[2])
out[2] = qout[2]
print "%g + %g*x^%g"%(out[0],out[1],out[2])
产量
40.1253 + -282.949 * x ^ 0.375555
关于python - 试图从scipy powerlaw fit中获得合理的值,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/10181151/