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 よりも大きな分数の電力依存性を見つけることを期待していました。

私は自分のデータを再スケーリングし、optimize.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 に答える