试图从scipy powerlaw fi获得合理的价值

Question

我试着从我运行的模拟代码中获取一些数据，以找出幂律依赖性。当我绘制线性拟合图时，数据拟合得不太好。

下面是我用来调整数据的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

虽然在情节上的配合看起来和你期望的一样好，从一个租赁方配合，输出的形式困扰我。我希望常数接近你期望的0（大约30）。我希望找到一个大于10^5的幂依赖。

我试着重新调整我的数据，玩参数优化.leatsq没有运气。我要完成的是可能的还是我的数据不允许？计算是昂贵的，因此获得更多的数据点是非常重要的。

谢谢！

Answer 1

使用线性最小二乘法获得指数拟合的标准方法是做fraxel suggests in his/her answer：将直线拟合到对数（y）。

然而，这种方法有已知的数值缺点，特别是灵敏度（数据的小变化会导致估计值的大变化）。首选的方法是使用非线性最小二乘法——它不太敏感。但如果您对线性LS方法的非关键用途感到满意，只需使用它。

Answer 2

它有助于重新缩放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

Answer 3

最好先取对数，然后用leastsquare来拟合这个线性方程，这样可以得到更好的拟合。在scipy cookbook中有一个很好的例子，我在下面对其进行了修改，以适合您的代码。

这样的最佳拟合是：振幅=0.8955，指数=0.40943265484

从图（和你的数据）中我们可以看到，如果它符合幂律，我们就不会期望振幅值接近30。如幂律方程f(x) == Amp * x ** index所示，负指数：f(1) == Amp和f(0) == infinity。

from pylab import *
from scipy import *
from scipy import optimize

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]

logx = log10(xdata)
logy = log10(ydata)

# define our (line) fitting function
fitfunc = lambda p, x: p[0] + p[1] * x
errfunc = lambda p, x, y: (y - fitfunc(p, x))

pinit = [1.0, -1.0]
out = optimize.leastsq(errfunc, pinit,
                       args=(logx, logy), full_output=1)

pfinal = out[0]
covar = out[1]

index = pfinal[1]
amp = 10.0**pfinal[0]

print 'amp:',amp, 'index', index

powerlaw = lambda x, amp, index: amp * (x**index)
##########
# Plotting data
##########
clf()
subplot(2, 1, 1)
plot(xdata, powerlaw(xdata, amp, index))     # Fit
plot(xdata, ydata)#, yerr=yerr, fmt='k.')  # Data
text(0.0020, 30, 'Ampli = %5.2f' % amp)
text(0.0020, 25, 'Index = %5.2f' % index)
xlabel('X')
ylabel('Y')

subplot(2, 1, 2)
loglog(xdata, powerlaw(xdata, amp, index))
plot(xdata, ydata)#, yerr=yerr, fmt='k.')  # Data
xlabel('X (log scale)')
ylabel('Y (log scale)')

savefig('power_law_fit.png')
show()