Fitting Lower Incomplete Gamma Models in Python: Avoid SymPy lambdify Pitfalls with SciPy and lmfit, plus a Clean Parameter Constraint

Learn to fit lower incomplete gamma models in Python with lmfit and SciPy, avoid SymPy lambdify errors, and implement parameter constraints as derived values.

Fitting special-function models in Python often collides with the boundaries between symbolic math and numeric evaluation. A frequent trap is mixing SymPy expressions with NumPy-backed lambdify calls when the target function is not properly supported. Below is a focused walkthrough on fitting a model that includes the lower incomplete gamma function with lmfit, plus a parameter constraint derived from the equation ((A*m**k)/k)-B=1.

Problem setup

The target model is a sum of two terms: a scaled lower incomplete gamma component and an exponential term. The goal is to fit this model to data using lmfit and to enforce a constraint tying parameters together. The original attempt uses SymPy to build the expressions and lambdify to generate a numeric function for lmfit.

Code that reproduces the issue

The example below mirrors the structure of the problematic approach. It builds the model with SymPy, lambdifies it, and sets up lmfit with a constraint. Note the naming and structure are consistent in logic with the described problem.

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from lmfit import Model, Parameters
from sympy.parsing import sympy_parser as sp_parser
# Data
x_points = [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]
y_points = [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]
# Symbolic model using lowergamma
sym_gamma = sp_parser.parse_expr('C*lowergamma(r,u*x)/(x**r)')
sym_exp = sp_parser.parse_expr('-1*D*exp(-v*x)')
parts = sp.Array((sym_gamma, sym_exp))
sym_expr = sum(parts)
# Lambdify for numerical evaluation
parts_fn = sp.lambdify(list(parts.free_symbols), parts)
expr_fn = sp.lambdify(list(sym_expr.free_symbols), sym_expr)
# lmfit model
generic_model = Model(expr_fn)
print(f'parameter names: {generic_model.param_names}')
print(f'independent variables: {generic_model.independent_vars}')
# Parameters and an additional expression intended as a constraint
fit_params = generic_model.make_params(D=dict(value=1, min=0, max=15),
                                       r=dict(value=7.2E-1, min=0),
                                       C=dict(value=.5, min=0, max=50),
                                       v=dict(value=1.31E-2, min=0),
                                       u=dict(value=1.31E-2, min=0))
fit_params.add('constraint_1', expr='((C*u**r)/r)-D', min=.9, max=1.1)
# Attempt to fit; note that PLt is referenced
result_bad = generic_model.fit(PLt[1], params=fit_params, x=PLt[0])

What goes wrong and why

The core blocker is that sympy.lowergamma is not recognized properly in the NumPy domain used by lambdify. This prevents the lambdified function from evaluating the lower incomplete gamma term as expected. The downstream fit cannot proceed reliably because the generated numeric function is not fully supported.

There is also a practical inconsistency in the fitting call: the code passes PLt, which is not defined in the snippet, instead of the provided data arrays. This invites additional errors unrelated to the special-function issue.

Finally, the parameter relation ((A*m**k)/k)-B=1 implies B is not an independent degree of freedom. Enforcing it as a separate bounded parameter is unnecessary if B is deterministically derived from A, m, and k.

Working solution with the same model

To make the model numerically evaluable, use SciPy’s special functions. The lower incomplete gamma can be expressed via gammainc and gamma from scipy.special. To honor the parameter relation, define B as a derived quantity inside the model function instead of fitting B explicitly.

from lmfit import Model, Parameters, report_fit
import numpy as np
from scipy.special import gammainc, gamma as gamma_fn
import matplotlib.pyplot as plt
# Data
t_vals = 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])
obs_vals = 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])
# Lower incomplete gamma via SciPy
def inc_gamma(shape, arg):
    return gammainc(shape, arg) * gamma_fn(shape)
# Model with B derived from the relation ((A*m**k)/k)-B=1
# A, m, k, q are the free parameters
# B = (A * m**k) / k - 1
def composite_model(t, alpha, kappa, mu, theta):
    beta = (alpha * mu**kappa) / kappa - 1
    g_part = alpha * inc_gamma(kappa, mu * t) / (t**kappa)
    e_part = -beta * np.exp(-theta * t)
    return g_part + e_part
# Build and fit
fit_obj = Model(composite_model)
init = Parameters()
init.add('alpha', value=1.0, min=0.01, max=10)
init.add('mu',    value=0.01, min=1e-5, max=1.0)
init.add('theta', value=0.01, min=1e-5, max=1.0)
init.add('kappa', value=0.72, min=0.1,  max=10)
fit_ok = fit_obj.fit(obs_vals, init, t=t_vals)
report_fit(fit_ok)
# Visualization
plt.plot(t_vals, obs_vals, 'o', label='data')
plt.plot(t_vals, fit_ok.best_fit, '-', label='fit')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.title('Fitted Model')
plt.grid(True)
plt.show()

Why this matters

When a special function is not available in the numeric backend used by lambdify, everything downstream becomes fragile. Replacing the symbolic special function with a SciPy implementation aligns the entire pipeline with NumPy’s array semantics and ensures lmfit can evaluate the model during optimization. Treating the parameter relation as a true constraint by deriving B from the other parameters also improves identifiability and avoids redundant degrees of freedom.

Practical notes and takeaways

If you pass arrays to lmfit, make sure the variables you reference in the fit call match the data you prepared; a stray name like PLt will derail the run immediately. When debugging, always include the full error message; it usually carries the vital hint pointing to the unsupported function or a mismatched argument. If the resulting percent error seems high, consider that it may stem from the chosen functional form and the number of parameters rather than from the mechanics of the fit itself.

The bottom line is straightforward. Keep special functions numeric by using SciPy where possible, eliminate redundant parameters by turning algebraic constraints into definitions inside the model function, and verify that your data variables match what you pass to the optimizer.