使用lmfit拟合含下不完全伽马函数模型时遇NaN错误及参数约束求助

尝试用lmfit将包含下不完全伽马函数的自定义函数拟合到给定数据集，同时施加参数约束 ((A*m**k)/k)-B=1，但反复触发错误：

ValueError: The model function generated NaN values and the fit aborted!

用户提供的代码如下：

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

问题分析及修复方案 #

1. 参数约束实现错误 #

当前添加约束的方式无法真正限制参数关系，反而会引入无效的派生参数。正确做法是通过表达式替代参数，减少自由参数：
根据约束 ((A*m**k)/k)-B=1，解出 B = (A*m**k)/k - 1，将B设为派生参数而非自由参数：

params = lm_mod.make_params(k=dict(value=7.2E-1, min=1e-6, max=10), 
                            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('B', expr='(A*m**k)/k - 1')

2. k参数的下限问题 #

原代码中k的下限设为0，当拟合过程中k趋近于0时，lowergamma(k, m*t)和分母k会引发NaN。给k设置极小的非零下限：

k=dict(value=7.2E-1, min=1e-6, max=10)

3. 数据传入错误 #

代码最后一行使用了未定义的PLt，应替换为实际定义的xdata和ydata：

res = lm_mod.fit(ydata, params=params, t=xdata)

4. 数值计算稳定性优化 #

SymPy的lowergamma转换为数值函数时稳定性不足，建议直接使用SciPy的专用函数实现模型，避免解析过程的潜在问题：

from scipy.special import gammainc, gamma

def model_func(t, A, k, m, B, q):
    # lowergamma(k, x) = gamma(k)*gammainc(k, x)
    term1 = A * gamma(k) * gammainc(k, m*t) / (t**k)
    term2 = -B * np.exp(-q*t)
    return term1 + term2

lm_mod = Model(model_func)

5. t**k的数值溢出问题 #

当t很小且k为极大值时，t**k可能下溢为0，引发NaN。可以将其改写为指数形式避免：

term1 = A * gamma(k) * gammainc(k, m*t) / np.exp(k * np.log(t))

内容的提问来源于stack exchange，提问作者MOC