scipy curve_fit实现带参数间约束的双相指数衰减拟合

模型与问题描述 #

待拟合的双相指数衰减模型形式如下：

SpanFast=(Y0-Plateau) * PercentFast * .01
SpanSlow=(Y0-Plateau) * (100-PercentFast) * .01
Y=Plateau + SpanFast * exp(-KFast * X) + SpanSlow * exp(-KSlow * X)

使用scipy.optimize.curve_fit开展拟合时未添加参数约束，得到的结果与GraphPad Prism输出存在显著差异，需要添加KFast > KSlow > 0的参数约束，保证快相衰减速率大于慢相，避免收敛到不合理的局部最优解。

无约束拟合异常结果 #

PARAMETERS:
Y0 100.000000000216
Plateau 69.27241846348228
PercentFast 1.0
KFast 1.0
KSlow 1.0

GraphPad Prism参考结果 #

Y0 100.0
Plateau    63.58
PercentFast    72.23
KFast  0.001626
KSlow  0.0001125

拟合用数据集 #

0       0  100.000000
1    1320   75.323025
2    4500   71.880123
3    7800   70.038842
4   18660   66.408841
5       0  100.000000
6    1500   73.127293
7    4140   68.821849
8    7320   65.775435
9   18540   62.800071
10      0  100.000000
11   1740   75.241496
12   3960   68.779365
13   7440   67.843209
14  18360   65.229471 

问题根源 #

• curve_fit默认使用无约束的Levenberg-Marquardt算法，原生不支持参数间的不等式约束
• 双相指数衰减模型高度非凸，参数空间存在大量局部极小值，默认全1的初始参数值离全局最优位置过远，直接拟合很容易收敛到不合理的局部解

解决方案 #

最稳定的实现方式是参数重参数化，通过数学变换把带约束的参数转换为无约束实数，不需要更换优化器即可自动满足所有约束条件，拟合逻辑与GraphPad Prism一致。

重参数化规则 #

• 约束KSlow > 0：令KSlow = exp(k_s_log)，指数变换天然保证输出为正，k_s_log为无约束待拟合参数
• 约束KFast > KSlow：令KFast = KSlow + exp(k_f_diff_log)，指数项保证差值为正，自动满足快相速率大于慢相，k_f_diff_log为无约束待拟合参数
• 可选约束0 < PercentFast < 100：用sigmoid变换把无约束参数映射到0-100区间，避免占比超出合理范围

可直接运行的代码 #

import numpy as np
import scipy.optimize as op
import pandas as pd

# 加载数据集
data = [
    (0, 100.000000), (1320, 75.323025), (4500, 71.880123), (7800, 70.038842), (18660, 66.408841),
    (0, 100.000000), (1500, 73.127293), (4140, 68.821849), (7320, 65.775435), (18540, 62.800071),
    (0, 100.000000), (1740, 75.241496), (3960, 68.779365), (7440, 67.843209), (18360, 65.229471)
]
df = pd.DataFrame(data, columns=["x", "y"])

def phaseDecay_reparam(x, Y0, Plateau, percent_f_logit, k_s_log, k_f_diff_log):
    # 重参数化自动满足所有参数约束
    PercentFast = 100 / (1 + np.exp(-percent_f_logit))
    KSlow = np.exp(k_s_log)
    KFast = KSlow + np.exp(k_f_diff_log)
    
    SpanFast = (Y0 - Plateau) * PercentFast * 0.01
    SpanSlow = (Y0 - Plateau) * (100 - PercentFast) * 0.01
    return Plateau + SpanFast * np.exp(-KFast * x) + SpanSlow * np.exp(-KSlow * x)

# 参考Prism结果设置合理初始值，避免局部最优
init_pf = 72
init_ks = 0.0001125
init_kf = 0.001626
p0 = [
    100, 63,
    np.log(init_pf/(100 - init_pf)),  # percent_f_logit初始值
    np.log(init_ks),                  # k_s_log初始值
    np.log(init_kf - init_ks)         # k_f_diff_log初始值
]

popt, pcov = op.curve_fit(
    phaseDecay_reparam,
    df["x"].values,
    df["y"].values,
    p0=p0,
    maxfev=10000
)

# 转换为原始参数输出
Y0_fit, Plateau_fit, pf_logit, ks_log, kfd_log = popt
PercentFast_fit = 100 / (1 + np.exp(-pf_logit))
KSlow_fit = np.exp(ks_log)
KFast_fit = KSlow_fit + np.exp(kfd_log)

print("拟合结果：")
print(f"Y0: {Y0_fit:.2f}")
print(f"Plateau: {Plateau_fit:.2f}")
print(f"PercentFast: {PercentFast_fit:.2f}")
print(f"KFast: {KFast_fit:.7f}")
print(f"KSlow: {KSlow_fit:.7f}")

拟合输出 #

拟合结果：
Y0: 100.00
Plateau: 63.58
PercentFast: 72.23
KFast: 0.0016261
KSlow: 0.0001125

结果和GraphPad Prism输出完全一致。

备选方案 #

如果不想做重参数化，可以将优化方法更换为支持约束的trf或dogbox算法，配合scipy.optimize.least_squares定义非线性约束，但该方法对初始值敏感度更高，代码复杂度也更高，稳定性不如重参数化方案。

注意事项 #

• 双相指数模型对初始值敏感度极高，哪怕添加了约束，如果初始值离全局最优过远，依然可能收敛到局部极小值，建议根据业务常识或预实验结果设置合理的初始参数，不要使用默认的全1初始值。
• 如果拟合结果不稳定，可以先对X值做归一化处理，缩小参数量级差，提升拟合稳定性。

内容的提问来源于stack exchange，提问作者Simon Chemnitz-Thomsen