Python中拟合双指数函数时约束y0+a1+a2=1的实现方法

问题背景 #

需要将数据拟合到公式 y(t) = y0 + a1*e^(-k1*t) + a2*e^(-k2*t)，且满足约束：

• y0、a1、a2 均为正值
• y0 + a1 + a2 = 1

原代码尝试在拟合函数内赋值 y0=1-a1-a2，但因仍将y0作为独立拟合参数，导致约束未生效。

解决思路 #

核心是减少拟合参数数量，直接用 y0 = 1 - a1 - a2 替代独立参数y0，让拟合过程只优化a1、k1、a2、k2四个参数，同时通过参数边界保证：

• a1 > 0，a2 > 0
• 1 - a1 - a2 > 0（通过设置a1和a2的上界为1，配合合理的初始值实现）

修改后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

nai = np.array([
    [1.22212386e+00, 3.29623401e+01, 6.47035556e+01, 9.64437723e+01,
     1.28183988e+02, 1.59925204e+02, 1.91665421e+02, 2.23405636e+02,
     2.55145451e+02, 2.86886261e+02, 3.18627075e+02, 3.50366878e+02,
     3.82106686e+02, 4.13846500e+02, 4.45586303e+02, 4.77327118e+02,
     5.09067921e+02, 5.40807729e+02, 5.72547538e+02, 6.04287346e+02,
     6.36027155e+02, 6.67767964e+02, 6.99507773e+02, 7.31248581e+02,
     7.62988390e+02, 7.94728199e+02, 8.26468997e+02, 8.58208596e+02,
     8.89948201e+02, 9.21687811e+02, 9.53428532e+02, 9.85169347e+02,
     1.01690915e+03, 1.04864896e+03, 1.08038877e+03, 1.11212858e+03,
     1.14386839e+03, 1.17561019e+03, 1.20735000e+03, 1.23908982e+03,
     1.27082962e+03, 1.30256943e+03, 1.33431024e+03, 1.36605004e+03,
     1.39779085e+03, 1.42953066e+03, 1.46127048e+03, 1.49301128e+03,
     1.52475109e+03, 1.55649076e+03, 1.58823037e+03, 1.61997097e+03,
     1.65171057e+03, 1.68345118e+03, 1.71519078e+03, 1.74693039e+03,
     1.77866999e+03, 1.81040960e+03],
    [1., 0.92072946, 0.93392534, 0.92599632, 0.84650909, 0.89966649,
     0.85375192, 0.8359797, 0.81487989, 0.83947038, 0.78711305, 0.77231505,
     0.74861849, 0.67860406, 0.73989558, 0.78446446, 0.72500833, 0.72903779,
     0.69205116, 0.73345758, 0.71143277, 0.71248356, 0.67470885, 0.67626045,
     0.63480034, 0.68824323, 0.69961798, 0.68302773, 0.64354886, 0.64451215,
     0.63644302, 0.63616151, 0.63049777, 0.62942254, 0.67511705, 0.62839769,
     0.61536734, 0.58931514, 0.60952486, 0.59674224, 0.55770579, 0.58915328,
     0.5567563, 0.58849712, 0.50377624, 0.50890445, 0.5446239, 0.51436443,
     0.52945348, 0.55326716, 0.51701847, 0.5240366, 0.54487462, 0.50567212,
     0.49233296, 0.5107218, 0.56567484, 0.48162376]
]).astype(np.float64)

# 定义带约束的双指数衰减函数：y0 = 1 - a1 - a2，不将y0作为拟合参数
def biexponential_decay(t, a1, k1, a2, k2):
    y0 = 1 - a1 - a2
    return a1 * np.exp(-k1 * t) + a2 * np.exp(-k2 * t) + y0

# 设置参数边界：a1>0, k1>0, a2>0, k2>0；同时a1<1, a2<1（保证y0=1-a1-a2>0）
bounds = ([0, 0, 0, 0], [1, np.inf, 1, np.inf])

# 初始参数选择：a1+a2<1，确保初始y0为正
p0 = (0.5, 0.001, 0.3, 0.001)

# 执行拟合
params, covariance = curve_fit(biexponential_decay, nai[0], nai[1], p0=p0, bounds=bounds, method='trf')

# 提取拟合参数并计算y0
a1_fit, k1_fit, a2_fit, k2_fit = params
y0_fit = 1 - a1_fit - a2_fit

# 计算残差、MSE和R²
residuals = nai[1] - biexponential_decay(nai[0], a1_fit, k1_fit, a2_fit, k2_fit)
mse = np.mean(residuals**2)
r_squared = 1 - mse / np.var(nai[1])

# 计算标准误差
cov_diag = np.diag(covariance)
std_error = np.sqrt(mse * cov_diag)

# 输出结果
print("拟合参数：")
print(f"a1 = {a1_fit:.6f}")
print(f"k1 = {k1_fit:.6f}")
print(f"a2 = {a2_fit:.6f}")
print(f"k2 = {k2_fit:.6f}")
print(f"y0 = {y0_fit:.6f}")
print(f"约束验证：y0+a1+a2 = {y0_fit+a1_fit+a2_fit:.6f}")
print(f"R² = {r_squared:.6f}")
print(f"参数标准误差：{std_error}")

# 绘制拟合结果
y_fit = biexponential_decay(nai[0], a1_fit, k1_fit, a2_fit, k2_fit)
plt.scatter(nai[0], nai[1], label='原始数据')
plt.plot(nai[0], y_fit, 'r-', label='拟合曲线')
plt.legend()
plt.xlabel('t')
plt.ylabel('y(t)')
plt.show()

关键修改说明 #

1. 函数参数简化：移除了y0作为拟合参数，直接在函数内通过y0=1-a1-a2计算，从根源上保证和为1的约束。
2. 边界调整：
   • a1和a2的下界设为0，上界设为1，确保两者为正且1-a1-a2也为正
   • k1和k2的下界设为0（衰减常数应为正），上界设为无穷大，保留足够的拟合空间
3. 初始值设置：选择a1=0.5、a2=0.3，初始y0=0.2为正，帮助拟合算法快速收敛到符合约束的解。

效果验证 #

运行代码后，输出的约束验证项会显示y0+a1+a2=1.0，严格满足和为1的要求；同时a1、a2、y0均为正值，符合物理系统的约束条件。

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