11

我正在尝试从我一直在运行的模拟代码中拟合一些数据，以找出幂律依赖性。当我绘制线性拟合时，数据拟合得不是很好。

这是我用来拟合数据的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 更大的功率依赖性。

我尝试重新调整我的数据并使用参数来优化.leastsq，但没有运气。我想要完成的事情是可能的，还是我的数据不允许这样做？计算成本很高，因此获得更多数据点并非易事。

谢谢！

4

3 回答 3

8

最好先取对数，然后用它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()
于 2012-04-16T22:19:20.487 回答
4

它有助于重新调整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

于 2012-04-16T22:00:03.270 回答
1

使用线性最小二乘法获得指数拟合的标准方法是按照fraxel 在他/她的回答中建议的方法：将直线拟合到 log(y_i)。

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

于 2012-04-16T22:34:36.830 回答