3

多分誰かが私を助けることができます。何日も費やしましたが、問題を解決できませんでした。前もって感謝します。

実験データに 2 つのローレンツを当てはめたいと思います。方程式を 2 つのローレンツ行列lorentz1とlorentz2関数の単純な形に分解しました。次に、他の 2 つの関数を定義L1しL2、定数を乗算するだけcnstにしました。適合する 4 つのパラメータがすべてあります: cnst1、cnst2、tau1、tau2。

私が使用するものlmfit: モデルと最小化 (おそらく両方とも同じ方法を使用します)。

初期の適合パラメータは、視覚的に微適合に近づくように設定されています。ただし、lmfit を使用した最小化は失われます (下の最初の画像)。

制約のないパラメーター

これらのパラメータを使用して:

params.add('cnst1', value=1e3 , min=1e2, max=1e5)
params.add('cnst2', value=3e5, min=1e2, max=1e6)
params.add('tau1', value=2e0, min=0, max=1e2)
params.add('tau2', value=5e-3, min=0, max=10)

しかし、エラー率は低いです:

cnst1:   117.459806 +/- 14.67188 (12.49%) (init= 1000)
cnst2:   413.959032 +/- 44.21042 (10.68%) (init= 300000)
tau1:    11.0343531 +/- 1.065570 (9.66%) (init= 2)
tau2:    1.55259664 +/- 0.125853 (8.11%) (init= 0.005)

一方、パラメーターを初期値に非常に近づける (強制的に初期値に近づける):

ここに画像の説明を入力

パラメータの使用:

#params.add('cnst1', value=1e3 , min=0.1e3, max=1e3)
#params.add('cnst2', value=3e5, min=1e3, max=1e6)
#params.add('tau1', value=2e0, min=0, max=2)
#params.add('tau2', value=5e-3, min=0, max=10)

フィットは視覚的に優れていますが、エラー値は非常に大きくなります:

[[Variables]]
cnst1:   752.988629 +/- 221.3098 (29.39%) (init= 1000)
cnst2:   3.0159e+05 +/- 3.05e+07 (10113.40%) (init= 300000)
tau1:    1.99684317 +/- 0.600748 (30.08%) (init= 2)
tau2:    0.00497806 +/- 0.289651 (5818.56%) (init= 0.005)

総コードは次のとおりです。

import numpy as np
from lmfit import Model, minimize, Parameters, report_fit
import matplotlib.pyplot as plt

x = np.array([0.02988, 0.07057,0.19365,0.4137,0.91078,1.85075,3.44353,6.39428,\
        11.99302,24.37024,52.58804,121.71927,221.53799,358.27392,464.70405])

y = 1.0 / np.array([4.60362E-4,5.63559E-4,8.44538E-4,0.00138,0.00287,0.00657,0.01506,\
            0.03119,0.0584,0.09153,0.12538,0.19389,0.34391,0.68869,1.0])

def lorentz1(x, tau):
    L =  tau  / ( 1 + (x*tau)**2 )  
    return(L)

def lorentz2(x, tau):
    L =  tau**2  / ( 1 + (x*tau)**2 )  
    return(L)

def L1(x,cnst1,tau1):
    L1 =  cnst1 * lorentz1(x,tau1)
    return (L1)

def L2(x, cnst2, tau2):
    L2 =  cnst2 * lorentz2(x,tau2)
    return (L2)    

def L_min(params, x, y):
    cnst1 = params['cnst1'].value
    cnst2 = params['cnst2'].value
    tau1 = params['tau1'].value
    tau2 = params['tau2'].value

    L_total = L1(x, cnst1, tau1) + L2(x, cnst2, tau2)
    resids = L_total - y
    return resids

#params  = mod.make_params( cnst1=10e2, cnst2=3e5, tau1=2e0, tau2=0.5e-2)
params = Parameters()
#params.add('cnst1', value=1e3 , min=0.1e3, max=1e3)
#params.add('cnst2', value=3e5, min=1e3, max=1e6)
#params.add('tau1', value=2e0, min=0, max=2)
#params.add('tau2', value=5e-3, min=0, max=10)

params.add('cnst1', value=1e3 , min=1e2, max=1e5)
params.add('cnst2', value=3e5, min=1e2, max=1e6)
params.add('tau1', value=2e0, min=0, max=1e2)
params.add('tau2', value=5e-3, min=0, max=10)


#1-----Model--------------------
mod = Model(L1) + Model(L2)
result_mod = mod.fit(y, params, x=x)
print('---results from lmfit.Model----')
print(result_mod.fit_report())

#2---minimize-----------
result_min = minimize(L_min, params, args=(x,y))
final_min = y + result_min.residual
print('---results from lmfit.minimize----')
report_fit(params)

#-------Plot------
plt.close('all')
plt.loglog(x, y,'bo' , label='experimental data')
plt.loglog(x, result_mod.init_fit, 'k--', label='initial')
plt.loglog(x, result_mod.best_fit, 'r-', label='final')
plt.legend()
plt.show()
4

1 に答える 1

4

何かを検索していると、Google は私が以前に尋ねた自分自身の質問をもたらしました。今、私は答えを知っているので、ここでそれを示します。誰かの役に立てば幸いです。:)

機能を考えてみlmfit.minimizeます。したがって、私が行った変更は、 の結果をプロットすることでしlmfit.minimizeた。そして、対数 y スケールの問題 (これは @mdurant によっても言及された主な問題でした) を解決するために、残差をその y 値に分割しました (何らかの方法で、残差を取得するときにすべてのデータを比較できるように正規化します)。私はこれを加重残差と名付けました。

def L_min(params, x, y):
    ...
    ..
    .
    resids = L_total - y
    weighted_resids = resids/y
    return weighted_resids

したがって、結果は青い線で示されます。

新しいフィッティング

そして完全なコード:

import numpy as np
from lmfit import Model, minimize, Parameters, report_fit
import matplotlib.pyplot as plt

x = np.array([0.02988, 0.07057,0.19365,0.4137,0.91078,1.85075,3.44353,6.39428,\
        11.99302,24.37024,52.58804,121.71927,221.53799,358.27392,464.70405])

y = 1.0 / np.array([4.60362E-4,5.63559E-4,8.44538E-4,0.00138,0.00287,0.00657,0.01506,\
            0.03119,0.0584,0.09153,0.12538,0.19389,0.34391,0.68869,1.0])

def lorentz1(x, tau):
    L =  tau  / ( 1 + (x*tau)**2 )  
    return(L)

def lorentz2(x, tau):
    L =  tau**2  / ( 1 + (x*tau)**2 )  
    return(L)

def L1(x,cnst1,tau1):
    L1 =  cnst1 * lorentz1(x,tau1)
    return (L1)

def L2(x, cnst2, tau2):
    L2 =  cnst2 * lorentz2(x,tau2)
    return (L2)    

def L_min(params, x, y):
    cnst1 = params['cnst1'].value
    cnst2 = params['cnst2'].value
    tau1 = params['tau1'].value
    tau2 = params['tau2'].value

    L_total = L1(x, cnst1, tau1) + L2(x, cnst2, tau2)
    resids = L_total - y
    weighted_resids = resids/y
    return weighted_resids
#    return resids

#params  = mod.make_params( cnst1=10e2, cnst2=3e5, tau1=2e0, tau2=0.5e-2)
params = Parameters()
#params.add('cnst1', value=1e3 , min=0.1e3, max=1e3)
#params.add('cnst2', value=3e5, min=1e3, max=1e6)
#params.add('tau1', value=2e0, min=0, max=2)
#params.add('tau2', value=5e-3, min=0, max=10)

params.add('cnst1', value=1e3 , min=1e2, max=1e5)
params.add('cnst2', value=3e5, min=1e2, max=1e6)
params.add('tau1', value=2e0, min=0, max=1e2)
params.add('tau2', value=5e-3, min=0, max=10)


#1-----Model--------------------
mod = Model(L1) + Model(L2)
result_mod = mod.fit(y, params, x=x)
print('---results from lmfit.Model----')
print(result_mod.fit_report())

#2---minimize-----------
result_min = minimize(L_min, params, args=(x,y))
final_min = y + result_min.residual
print('---results from lmfit.minimize----')
report_fit(params)

#-------Plot------
plt.close('all')
plt.loglog(x, y,'bo' , label='experimental data')
plt.loglog(x, result_mod.init_fit, 'k--', label='initial')
plt.loglog(x, result_mod.best_fit, 'r-', label='lmfit.Model')

min_result = L1(x, params['cnst1'].value, params['tau1'].value) + \
            L2(x, params['cnst2'].value, params['tau2'].value)
plt.loglog(x, min_result, 'b-', label='lmfit.Minimize')
plt.legend()
plt.show()

フィッティングエラーは素晴らしいです：

    cnst1:   832.592441 +/- 77.32939 (9.29%) (init= 1000)
    cnst2:   2.0836e+05 +/- 3.55e+04 (17.04%) (init= 300000)
    tau1:    1.64355457 +/- 0.221466 (13.47%) (init= 2)
    tau2:    0.00700899 +/- 0.000935 (13.34%) (init= 0.005)
[[Correlations]] (unreported correlations are <  0.100)
    C(tau1, tau2)                =  0.151
于 2015-01-09T14:42:09.037 に答える