2016 年 4 月 6 日更新
ほとんどの場合、当てはめパラメーターで正しいエラーを取得することは微妙な場合があります。
y=f(x)一連のデータ ポイントがある関数の当てはめについて考えてみましょう(x_i, y_i, yerr_i)。ここで、iは各データ ポイントにまたがるインデックスです。
ほとんどの物理的測定では、エラーyerr_iは測定デバイスまたは手順の系統的な不確実性であるため、に依存しない定数と考えることができますi。
どのフィッティング関数を使用するか、およびパラメーター エラーを取得する方法は?
optimize.leastsqメソッドは分数共分散行列を返します。この行列のすべての要素に残差分散 (つまり、減らされたカイ 2 乗) を掛けて、対角要素の平方根を取ると、適合パラメーターの標準偏差が推定されます。そのためのコードを以下の関数の 1 つに含めました。
一方、 を使用する場合optimize.curvefitは、上記の手順の最初の部分 (減らされたカイ 2 乗による乗算) が舞台裏で実行されます。次に、共分散行列の対角要素の平方根をとって、適合パラメーターの標準偏差を推定する必要があります。
さらに、データ ポイントごとoptimize.curvefitに値が異なる、より一般的なケースに対処するためのオプションのパラメーターを提供します。ドキュメントyerr_iから：
sigma : None or M-length sequence, optional
If not None, the uncertainties in the ydata array. These are used as
weights in the least-squares problem
i.e. minimising ``np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )``
If None, the uncertainties are assumed to be 1.
absolute_sigma : bool, optional
If False, `sigma` denotes relative weights of the data points.
The returned covariance matrix `pcov` is based on *estimated*
errors in the data, and is not affected by the overall
magnitude of the values in `sigma`. Only the relative
magnitudes of the `sigma` values matter.
エラーが正しいことを確認するにはどうすればよいですか?
適合パラメータの標準誤差の適切な推定値を決定することは、複雑な統計問題です。共分散行列の結果は、エラーの確率分布とパラメーター間の相互作用に関する仮定によって実装されoptimize.curvefit、実際に依存しています。optimize.leastsq特定の適合関数に応じて、存在する可能性のある相互作用f(x)。
私の意見では、複雑な問題に対処する最善の方法は、このリンクf(x)で概説されているブートストラップ メソッドを使用することです。
いくつかの例を見てみましょう
まず、いくつかのボイラープレート コードです。波線関数を定義して、ランダム エラーを含むデータを生成してみましょう。小さなランダムエラーでデータセットを生成します。
import numpy as np
from scipy import optimize
import random
def f( x, p0, p1, p2):
return p0*x + 0.4*np.sin(p1*x) + p2
def ff(x, p):
return f(x, *p)
# These are the true parameters
p0 = 1.0
p1 = 40
p2 = 2.0
# These are initial guesses for fits:
pstart = [
p0 + random.random(),
p1 + 5.*random.random(),
p2 + random.random()
]
%matplotlib inline
import matplotlib.pyplot as plt
xvals = np.linspace(0., 1, 120)
yvals = f(xvals, p0, p1, p2)
# Generate data with a bit of randomness
# (the noise-less function that underlies the data is shown as a blue line)
xdata = np.array(xvals)
np.random.seed(42)
err_stdev = 0.2
yvals_err = np.random.normal(0., err_stdev, len(xdata))
ydata = f(xdata, p0, p1, p2) + yvals_err
plt.plot(xvals, yvals)
plt.plot(xdata, ydata, 'o', mfc='None')
それでは、利用可能なさまざまな方法を使用して関数を当てはめましょう。
`optimize.leastsq`
def fit_leastsq(p0, datax, datay, function):
errfunc = lambda p, x, y: function(x,p) - y
pfit, pcov, infodict, errmsg, success = \
optimize.leastsq(errfunc, p0, args=(datax, datay), \
full_output=1, epsfcn=0.0001)
if (len(datay) > len(p0)) and pcov is not None:
s_sq = (errfunc(pfit, datax, datay)**2).sum()/(len(datay)-len(p0))
pcov = pcov * s_sq
else:
pcov = np.inf
error = []
for i in range(len(pfit)):
try:
error.append(np.absolute(pcov[i][i])**0.5)
except:
error.append( 0.00 )
pfit_leastsq = pfit
perr_leastsq = np.array(error)
return pfit_leastsq, perr_leastsq
pfit, perr = fit_leastsq(pstart, xdata, ydata, ff)
print("\n# Fit parameters and parameter errors from lestsq method :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from lestsq method :
pfit = [ 1.04951642 39.98832634 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
`optimize.curve_fit`
def fit_curvefit(p0, datax, datay, function, yerr=err_stdev, **kwargs):
"""
Note: As per the current documentation (Scipy V1.1.0), sigma (yerr) must be:
None or M-length sequence or MxM array, optional
Therefore, replace:
err_stdev = 0.2
With:
err_stdev = [0.2 for item in xdata]
Or similar, to create an M-length sequence for this example.
"""
pfit, pcov = \
optimize.curve_fit(f,datax,datay,p0=p0,\
sigma=yerr, epsfcn=0.0001, **kwargs)
error = []
for i in range(len(pfit)):
try:
error.append(np.absolute(pcov[i][i])**0.5)
except:
error.append( 0.00 )
pfit_curvefit = pfit
perr_curvefit = np.array(error)
return pfit_curvefit, perr_curvefit
pfit, perr = fit_curvefit(pstart, xdata, ydata, ff)
print("\n# Fit parameters and parameter errors from curve_fit method :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from curve_fit method :
pfit = [ 1.04951642 39.98832634 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
`ブートストラップ`
def fit_bootstrap(p0, datax, datay, function, yerr_systematic=0.0):
errfunc = lambda p, x, y: function(x,p) - y
# Fit first time
pfit, perr = optimize.leastsq(errfunc, p0, args=(datax, datay), full_output=0)
# Get the stdev of the residuals
residuals = errfunc(pfit, datax, datay)
sigma_res = np.std(residuals)
sigma_err_total = np.sqrt(sigma_res**2 + yerr_systematic**2)
# 100 random data sets are generated and fitted
ps = []
for i in range(100):
randomDelta = np.random.normal(0., sigma_err_total, len(datay))
randomdataY = datay + randomDelta
randomfit, randomcov = \
optimize.leastsq(errfunc, p0, args=(datax, randomdataY),\
full_output=0)
ps.append(randomfit)
ps = np.array(ps)
mean_pfit = np.mean(ps,0)
# You can choose the confidence interval that you want for your
# parameter estimates:
Nsigma = 1. # 1sigma gets approximately the same as methods above
# 1sigma corresponds to 68.3% confidence interval
# 2sigma corresponds to 95.44% confidence interval
err_pfit = Nsigma * np.std(ps,0)
pfit_bootstrap = mean_pfit
perr_bootstrap = err_pfit
return pfit_bootstrap, perr_bootstrap
pfit, perr = fit_bootstrap(pstart, xdata, ydata, ff)
print("\n# Fit parameters and parameter errors from bootstrap method :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from bootstrap method :
pfit = [ 1.05058465 39.96530055 1.96074046]
perr = [ 0.06462981 0.1118803 0.03544364]
観察
すでに興味深いことがわかり始めており、3 つの方法すべてのパラメーターと誤差の推定値はほぼ一致しています。それはいいです！
ここで、データに他の不確実性があることをフィッティング関数に伝えたいとします。おそらく、 の値の 20 倍の追加誤差に寄与する系統的不確実性ですerr_stdev。これは多くのエラーです。実際、そのようなエラーのあるデータをシミュレートすると、次のようになります。
このレベルのノイズで適合パラメーターを回復できる見込みはまったくありません。
まず、leastsqこの新しい系統誤差情報を入力することさえできないことを理解しましょう。curve_fitエラーについて伝えるとどうなるか見てみましょう。
pfit, perr = fit_curvefit(pstart, xdata, ydata, ff, yerr=20*err_stdev)
print("\nFit parameters and parameter errors from curve_fit method (20x error) :")
print("pfit = ", pfit)
print("perr = ", perr)
Fit parameters and parameter errors from curve_fit method (20x error) :
pfit = [ 1.04951642 39.98832633 1.95947613]
perr = [ 0.0584024 0.10597135 0.03376631]
なに?? これは絶対に間違っているに違いない！
これで話は終わりでしたが、最近、オプションのパラメーターがcurve_fit追加されました。absolute_sigma
pfit, perr = fit_curvefit(pstart, xdata, ydata, ff, yerr=20*err_stdev, absolute_sigma=True)
print("\n# Fit parameters and parameter errors from curve_fit method (20x error, absolute_sigma) :")
print("pfit = ", pfit)
print("perr = ", perr)
# Fit parameters and parameter errors from curve_fit method (20x error, absolute_sigma) :
pfit = [ 1.04951642 39.98832633 1.95947613]
perr = [ 1.25570187 2.27847504 0.72600466]
それはいくぶん良くなりましたが、それでも少し怪しいです。 curve_fitパラメータに 10% のレベルの誤差があり、そのノイズの多い信号から近似を得ることができると考えていp1ます。何を言わなければならないか見てみましょうbootstrap：
pfit, perr = fit_bootstrap(pstart, xdata, ydata, ff, yerr_systematic=20.0)
print("\nFit parameters and parameter errors from bootstrap method (20x error):")
print("pfit = ", pfit)
print("perr = ", perr)
Fit parameters and parameter errors from bootstrap method (20x error):
pfit = [ 2.54029171e-02 3.84313695e+01 2.55729825e+00]
perr = [ 6.41602813 13.22283345 3.6629705 ]
ああ、それはおそらく、適合パラメーターの誤差をより適切に見積もっています。約 34% の不確実性 bootstrapで知っていると考えています。p1
概要
optimize.leastsq当てはめたパラメーターの誤差を推定する方法をoptimize.curvefit提供しますが、これらの方法を少し疑問視せずに使用することはできません。これbootstrapは総当たりを使用する統計的方法であり、私の意見では、解釈が難しい状況でよりうまく機能する傾向があります.
特定の問題を見て、試してみることを強くお勧めcurvefitしbootstrapます。それらが類似している場合、curvefit計算がはるかに安価になるため、おそらく使用する価値があります。それらが大幅に異なる場合、私のお金はbootstrap.