Global fitting example with symfit

Question

I am trying to perform global fitting with 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()

In this example, I expect the "kg" parameter in the "globalmodel" is optimized to 0.005. However, the value of "kg" is about 9.6e-3, which is too near to the initial value (10.0e-3). I think I do something stupid, but I cannot figure it out.

Any comments and suggestions are welcome!

EDITED

I added (a very ugly) while loop to get the best fit. I am not sure why it should be, but it seems to work.

Answer 1

It would appear that the bounds are causing the problem. I removed them in my test and then everything works fine. 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.

Please note, I changed your np.e to sf.exp because that is symbolic. My working code is below, identical to yours except for the change mentioned and running in 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()