使用 symfit 的全局拟合示例

Global fitting example with symfit

我正在尝试使用 symfit package, following symfit documentation 执行全局拟合。

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 循环以获得最佳匹配。我不确定为什么会这样,但它似乎有效。

问题似乎是边界引起的。我在测试中删除了它们,然后一切正常。 This is a known problem in symfit 0.3.3, a̶n̶d̶ ̶o̶n̶e̶ ̶I̶ ̶a̶l̶r̶e̶a̶d̶y̶ ̶f̶i̶x̶e̶d̶ ̶i̶n̶ ̶t̶h̶e̶ ̶[̶̶m̶a̶s̶t̶e̶r̶̶]̶[̶1̶]̶ ̶b̶r̶a̶n̶c̶h̶ ̶o̶n̶ ̶G̶i̶t̶h̶u̶b̶.̶ ̶ ̶ ̶I̶ ̶u̶p̶l̶o̶a̶d̶e̶d̶ ̶a̶ ̶n̶e̶w̶ ̶d̶e̶v̶ ̶v̶e̶r̶s̶i̶o̶n̶ ̶y̶o̶u̶ ̶c̶o̶u̶l̶d̶ ̶n̶o̶w̶ ̶i̶n̶s̶t̶a̶l̶l̶ ̶u̶s̶i̶n̶g̶ ̶̶p̶i̶p̶ ̶i̶n̶s̶t̶a̶l̶l̶ ̶s̶y̶m̶f̶i̶t̶=̶=̶0̶.̶3̶.̶3̶.̶d̶e̶v̶1̶5̶5̶ ̶-̶-̶u̶p̶g̶r̶a̶d̶e̶̶,̶ ̶u̶n̶t̶i̶l̶ ̶I̶ ̶o̶f̶f̶i̶c̶i̶a̶l̶l̶y̶ ̶r̶e̶l̶e̶a̶s̶e̶ ̶0̶.̶3̶.̶4̶ ̶(̶w̶h̶i̶c̶h̶ ̶w̶i̶l̶l̶ ̶b̶e̶ ̶i̶d̶e̶n̶t̶i̶c̶a̶l̶ ̶b̶u̶t̶ ̶w̶i̶t̶h̶ ̶e̶x̶t̶e̶n̶d̶e̶d̶ ̶d̶o̶c̶u̶m̶e̶n̶t̶a̶t̶i̶o̶n̶)̶, which has now been fixed in newer versions.

请注意,我将您的 np.e 更改为 sf.exp,因为这是象征性的。我的工作代码如下,除了提到的更改和 0.3.3.dev155.

中的 运行 外,与您的相同
import numpy as np
import symfit as sf
import matplotlib.pyplot as plt

# 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 * sf.exp(- kg * x_1),
    y_2: a_2 * sf.exp(- 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)

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()