I am trying to fit a custom function, which involves a set lower incomplete gamma function, to a set of data I have using lmfit. I additionally want to impose a constraint on my parameters given by the equation ((A*m**k)/k)-B=1. I keep getting an error, saying it doesn't work.
Please help I do not know what I am doing.
Here is my code including the data.
import numpy as np
import sympy as sy
import matplotlib.pyplot as plt
from lmfit import Model
from lmfit import Parameters
from sympy.parsing import sympy_parser
#%%
xdata = [0.01, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360, 390, 420, 450, 480, 510, 540, 570, 600]
ydata = [1, 1.081070134, 1.123434136, 1.163246076, 1.180650102, 1.198810838, 1.20990884, 1.218026926, 1.221569822, 1.228107416, 1.223498562, 1.232861926, 1.23392966, 1.222959988, 1.21195955, 1.195828866, 1.174534424, 1.153189058, 1.136328876, 1.121319582, 1.100934304, 1.07886202, 1.073743626, 1.053117992, 1.035234418, 1.016680182, 1.002735456, 0.993092576, 0.97627191, 0.971622838, 0.94703108]
#%%
gamma = sympy_parser.parse_expr('A*lowergamma(k,m*t)/(t**k)')
sf = sympy_parser.parse_expr('-1*B*exp(-q*t)')
model_list = sy.Array((gamma, sf))
model = sum(model_list)
model_list_func = sy.lambdify(list(model_list.free_symbols), model_list)
model_func = sy.lambdify(list(model.free_symbols), model)
lm_mod = Model(model_func)
print(f'parameter names: {lm_mod.param_names}')
print(f'independent variables: {lm_mod.independent_vars}')
params = lm_mod.make_params(B=dict(value=1,min=0,max=15), k=dict(value=7.2E-1,min=0),
A=dict(value=.5,min=0,max=50), q=dict(value=1.31E-2,min=0), m=dict(value=1.31E-2,min=0))
params.add('constraint_1',expr='((A*m**k)/k)-B',min=.9,max=1.1)
res = lm_mod.fit(PLt[1], params=params, t=PLt[0])
sympy.lowergamma is not recognized properly in the NumPy domain used by lambdify. Check:https://github.com/sympy/sympy/issues/15134
Try using scipy package for lowergamma function. For the contrains, since you have ((A*m**k)/k)-B=1, then you do not need parameter B:
from lmfit import Model, Parameters, report_fit
import numpy as np
from scipy.special import gammainc, gamma as gamma_func
import matplotlib.pyplot as plt
# Your data
xdata = np.array([0.01, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360, 390, 420, 450, 480, 510, 540, 570, 600])
ydata = np.array([1, 1.081070134, 1.123434136, 1.163246076, 1.180650102, 1.198810838, 1.20990884, 1.218026926, 1.221569822, 1.228107416, 1.223498562, 1.232861926, 1.23392966, 1.222959988, 1.21195955, 1.195828866, 1.174534424, 1.153189058, 1.136328876, 1.121319582, 1.100934304, 1.07886202, 1.073743626, 1.053117992, 1.035234418, 1.016680182, 1.002735456, 0.993092576, 0.97627191, 0.971622838, 0.94703108])
# Define custom gamma function
def lowergamma(k, x):
return gammainc(k, x) * gamma_func(k)
# Model function with constrained B
def model_func(t, A, k, m, q):
B = (A * m**k) / k - 1
gamma_part = A * lowergamma(k, m * t) / (t**k)
exp_part = -B * np.exp(-q * t)
return gamma_part + exp_part
# lmfit model
lm_mod = Model(model_func)
# Parameters (no B here, it's derived)
params = Parameters()
params.add('A', value=1.0, min=0.01, max=10)
params.add('m', value=0.01, min=1e-5, max=1.0)
params.add('q', value=0.01, min=1e-5, max=1.0)
params.add('k', value=0.72, min=0.1, max=10)
# Fit the model
result = lm_mod.fit(ydata, params, t=xdata)
# Show fit report
report_fit(result)
# Plot
plt.plot(xdata, ydata, 'o', label='Data')
plt.plot(xdata, result.best_fit, '-', label='Best Fit')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.title("Fitted Model")
plt.grid(True)
plt.show()