symfit documentationに従って、 symfitパッケージでグローバル フィッティングを実行しようとしています。
import numpy as np
import symfit as sf
import matplotlib.pyplot as plt
%matplotlib inline # for ipynb
# Generate example data
t = np.arange(0.0, 600.1, 30)
k = 0.005
C1_0, C2_0 = 1.0, 2.0
C1 = C1_0 * np.exp(-k*t)
C2 = C2_0 * np.exp(-k*t)
# Construct model
x_1, x_2, y_1, y_2 = sf.variables('x_1, x_2, y_1, y_2')
kg = sf.Parameter(value=0.01, min=0.0, max=0.1)
a_1, a_2 = sf.parameters('a_1, a_2')
globalmodel = sf.Model({
y_1: a_1 * np.e**(- kg * x_1),
y_2: a_2 * np.e**(- kg * x_2),
})
# Do fit
globalfit = sf.Fit(globalmodel, x_1=t, x_2=t, y_1=C1, y_2=C2)
globalfit_result = globalfit.execute()
print(globalfit_result)
### EDITED START
while globalfit_result.r_squared < 0.99:
kg = sf.Parameter(value=globalfit_result.params['kg'])
a_1 = sf.Parameter(value=globalfit_result.params['a_1'])
a_2 = sf.Parameter(value=globalfit_result.params['a_2'])
globalmodel = sf.Model({
y_1: a_1 * np.e**(- kg * x_1),
y_2: a_2 * np.e**(- kg * x_2),
})
globalfit = sf.Fit(globalmodel, x_1=t, x_2=t, y_1=C1, y_2=C2)
globalfit_result = globalfit.execute()
### EDITED END
y_r = globalmodel(x_1=t, x_2=t, **globalfit_result.params)
# Plot fit
plt.plot(t,C1,'ro')
plt.plot(t,C2,'b+')
plt.plot(t,y_r[0],'r-')
plt.plot(t,y_r[1],'b-')
plt.show()
この例では、"globalmodel" の "kg" パラメータが 0.005 に最適化されていることを期待しています。ただし、「kg」の値は約 9.6e-3 であり、初期値 (10.0e-3) に近すぎます。私は愚かなことをしていると思いますが、それを理解することはできません。
コメントや提案は大歓迎です!
編集済み
最適なフィットを得るために (非常に醜い) while ループを追加しました。なぜそうなのかはわかりませんが、うまくいくようです。