python - 试图从scipy powerlaw fit中获得合理的值

我试图从我一直在运行的仿真代码中拟合一些数据，以找出幂律的相关性。当我绘制线性拟合时，数据拟合得不是很好。

这是我用来拟合数据的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/