Q&A
前提・実現したいこと
以下のサイトにあるサンプルを参考に、重み付き最小二乗法を用いてパラメータ推定を行いたいと考えております。
https://scipython.com/book/chapter-8-scipy/examples/weighted-and-non-weighted-least-squares-fitting/
サンプル
import numpy as np from scipy.optimize import curve_fit import pylab x0, A, gamma = 12, 3, 5 n = 200 x = np.linspace(1, 20, n) yexact = A * gamma**2 / (gamma**2 + (x-x0)**2) # Add some noise with a sigma of 0.5 apart from a particularly noisy region # near x0 where sigma is 3 sigma = np.ones(n)*0.5 sigma[np.abs(x-x0+1)<1] = 3 noise = np.random.randn(n) * sigma y = yexact + noise def f(x, x0, A, gamma): """ The Lorentzian entered at x0 with amplitude A and HWHM gamma. """ return A *gamma**2 / (gamma**2 + (x-x0)**2) def rms(y, yfit): return np.sqrt(np.sum((y-yfit)**2)) # Unweighted fit p0 = 10, 4, 2 popt, pcov = curve_fit(f, x, y, p0) yfit = f(x, *popt) print('Unweighted fit parameters:', popt) print('Covariance matrix:'); print(pcov) print('rms error in fit:', rms(yexact, yfit)) print() # Weighted fit popt2, pcov2 = curve_fit(f, x, y, p0, sigma=sigma, absolute_sigma=True) yfit2 = f(x, *popt2) print('Weighted fit parameters:', popt2) print('Covariance matrix:'); print(pcov2) print('rms error in fit:', rms(yexact, yfit2)) pylab.plot(x, yexact, label='Exact') pylab.plot(x, y, 'o', label='Noisy data') pylab.plot(x, yfit, label='Unweighted fit') pylab.plot(x, yfit2, label='Weighted fit') pylab.ylim(-1,4) pylab.legend(loc='lower center') pylab.show()
作成コード
現在は比較用に作っているUnweighted fitting用のコードしかかけていません。
import numpy as np import matplotlib.pyplot as plt import math from scipy.optimize import curve_fit # 係数見つけたい関数 def func(x, a): return ((-(a+(b*x))) + ((a+((b*x)**2)) - (4*b*math.log(0.1)))**0.5) / (2*b) # 求めたい分布関数の元データ x = np.array([6.26379, 8.57417, 8.66527, 8.75069, 11.6708, 12.3487, 14.5032, 15.7422, 21.7646, 23.0518, 26.5069, 26.4035, 26.321, 23.0045, 19.2654, 17.9425, 14.5669, 13.513, 10.4902, 9.95136, 9.77395]) y = np.array([3.709910308, 3.300454417, 3.219869361, 2.879991517, 2.250120678, 2.24981186, 1.859931899, 1.839996231, 1.560029151, 1.360016958, 1.210037387, 1.527926405, 1.320005022, 1.340038138, 1.618120234, 1.410033737, 1.83006856, 1.849465938, 2.141939621, 2.219958336, 2.494675074]) plt.plot(x, y, 'bo', label='Experimental data') plt.legend() # 観測データをフィッティングして係数を求める popt, pcov = curve_fit(func, x, y) popt # 各曲線の描画 plt.plot(x, y, 'bo', label="Experimental data") plt.plot(x, func(x, *popt),'ro', label="Unweighted fitting:a={:.3f}".format(*popt)) plt.legend()
Weighted fittingに関する箇所
ここを現在作成しているのですが、以下のコードをかくと、sigmaが定義されていないとエラーが出ます。x, yともにノイズを使用するのではなく、実験データを使用しようとした場合、どのようなコードを追加すれば宜しいのでしょうか？
# 観測データをフィッティングして係数を求める # WLS popt2, pcov2 = curve_fit(func, x, y, sigma=sigma, absolute_sigma=True) popt2
エラーは
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-44-d451ecf14c51> in <module> 2 # WLS 3 ----> 4 popt2, pcov2 = curve_fit(func, x, y, sigma=sigma, absolute_sigma=True) 5 popt2 NameError: name 'sigma' is not defined
言葉足らずではありますが、お力添えいただければ幸いです。
何卒宜しくお願いいたします。
回答1件
0
ベストアンサー
curve_fitの重み付き最小二乗法におけるsigmaは、まさにその「重み」を表すものなので、質問者様にて定義ください。
参考: Yamamoto's Laboratory - SciPyフィッティング (fitting)
sigma: None or M-length sequence or MxM array, optional デフォルト: sigma=None
ydataの不確かさの設定です．残差を r = ydata - f(xdata, popt). と定義すると，sigma の解釈は次元に依存します．
1次元の sigma は，y データのエラーの標準偏差です．この場合，最適化関数 (フィッティングにより最小値にする) は chisq = sum((r/sigma)* 2) です．
2次元の sigma は，y データの誤差の共分散行列です，最適化関数は chisq = r.T @ inv(sigma) @ r です．
None（デフォルト）は1で埋められた1-dシグマと同等です．
今回の場合、標準偏差sigmaをx, yから求める関数 sigma = s(x, y) : 数列×数列→数列 として定義すれば、xとyの観測値をそのまま入れることで、算出可能です。
ちなみに、sigma_scaler = s1(x, y) : スカラー×スカラー→スカラー として定義して、numpy.vectorizeを使えば、楽に、数列×数列→数列 の求める関数を作ることができます。
投稿2020/11/21 13:47
編集2020/11/21 13:55
総合スコア3266
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2020/11/21 14:06
2020/11/22 00:50 編集