gpt4 book ai didi

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

转载 作者:太空狗 更新时间:2023-10-29 18:03:16 25 4
gpt4 key购买 nike

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

这是我用来拟合数据的 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 但没有成功。我正在努力实现的目标是可能的还是我的数据不允许这样做?计算成本很高,因此获取更多数据点并非易事。

谢谢!

最佳答案

最好先取对数,然后使用leastsquare 来拟合这个线性方程,这会给你一个更好的拟合。 scipy cookbook 中有一个很好的例子,我在下面对其进行了调整以适合您的代码。

这样的最佳拟合是:amplitude = 0.8955,index = -0.40943265484

正如我们从图表(和您的数据)中看到的那样,如果它符合幂律,我们预计振幅值不会接近 30。正如在幂律方程 f(x) == Amp * x ** index 中,具有负指数:f(1) == Amp f(0) == 无穷大

enter image description here

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()

关于python - 试图从 scipy powerlaw fit 中获得合理的值(value),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/10181151/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com