在 python 中拟合双指数时，如何约束多个拟合参数总和为 1？

Question

我有一些衰减数据，看起来可能是双指数的，所以我想将其拟合到方程 y(t) = y0 + a1e^(-k1t) + a2e^(-k2t ），同时拟合 y0、a1、k1、a2 和 k2 的参数。我可以得到一个不错的拟合，但根据物理系统，我正在分析 y0、a1 和 a2 的参数总和应等于 1，并且所有三个均为正值。我有办法做到这一点吗？

我尝试过这个，它非常适合，但不会强制 a1 + a2 + y0 = 1。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
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]])

#%%biexponential decay
# Define the biexponential decay function with a1 + a2 + y0 = 1


nai=nai.astype(np.float64)

def biexponential_decay(nai, a1, k1, a2, k2, y0):
    y0 = 1 - a1 - a2
    return a1 * np.exp(-k1 * nai[0]) + a2 * np.exp(-k2 * nai[0]) + y0

# Define bounds for parameters
bounds = ([0, 0, 0, 0, 0], [1, 1, 1, 1, 1])  # Lower and upper bounds for a1, k1, a2, k2, and y0

# Fit the curve to the data
params, covariance = curve_fit(biexponential_decay, nai, nai[1], p0=(0.5, .001, 0.3, .001, .2), bounds=bounds, method='trf')  # Providing initial parameters as p0

# Extract the fitted parameters
a1_fit, k1_fit, a2_fit, k2_fit, y0_fit = params

# Calculate the residuals
residuals = nai[1] - biexponential_decay(nai, a1_fit, k1_fit, a2_fit, k2_fit, y0_fit)

# Calculate the mean squared error
mse = np.mean(residuals**2)

# Calculate the diagonal elements of the covariance matrix
cov_diag = np.diag(covariance)

# Calculate the standard error on the fitted values
std_error = np.sqrt(mse * cov_diag)

# Display the fitted parameters, R-squared value, and standard errors
print("Fitted Parameters:")
print("a1 =", a1_fit)
print("k1 =", k1_fit)
print("a2 =", a2_fit)
print("k2 =", k2_fit)
print("y0 =", y0_fit)
print("R-squared:", 1 - mse / np.var(nai[1]))
print("Standard Errors:", std_error)

# Generate the fitted curve using the fitted parameters
y_fit = biexponential_decay(nai, a1_fit, k1_fit, a2_fit, k2_fit, y0_fit)

# Plot the original data and the fitted curve
plt.scatter(nai[0], nai[1], label='Original Data')
plt.plot(nai[0], y_fit, 'r-', label='Fitted Curve')
plt.legend()
plt.xlabel('nai[0]')
plt.ylabel('nai[1]')
plt.show()

Answer 1

首先我会检查你是否真的需要 y0 为正值。我认为这是为了保持大 x 的正值，但您的数据实际上并没有接近这个值。您会得到相当好的曲线拟合，

a2

略大于 1，

y0

略为负值。

但是，如果您确实需要强制一个值，那么执行它的一个方法就是在不满足约束时向您的拟合函数添加一个巨大的惩罚。

在这里，我按照 mkrieger 的建议删除

y0

作为参数，只需将其设置为

1-a1-a2

即可。您应用了惩罚的函数是

def biexponential_decay(nai, a1, k1, a2, k2):
    y0 = 1 - a1 - a2
    penalty = 0
    if y0 < 0 : penalty = 1000000
    return a1 * np.exp(-k1 * nai[0]) + a2 * np.exp(-k2 * nai[0]) + y0 + penalty

Fitted Parameters:
a1 = 0.15346342751410955
k1 = 0.00647572574481048
a2 = 0.8465352091104801
k2 = 0.0002945049705031798
y0 = 1.363375410345924e-06
R-squared: 0.9582009677718569
Standard Errors: [1.60315165e-03 8.42485320e-05 2.00306402e-02 1.06429003e-05]

仅供参考，完整代码（已修改以删除未使用的导入并删除

y0

作为参数）是

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

#%%biexponential decay
# Define the biexponential decay function with a1 + a2 + y0 = 1


nai=nai.astype(np.float64)

def biexponential_decay(nai, a1, k1, a2, k2):
    y0 = 1 - a1 - a2
    penalty = 0
    if y0 < 0 : penalty = 1000000
    return a1 * np.exp(-k1 * nai[0]) + a2 * np.exp(-k2 * nai[0]) + y0 + penalty

# Define bounds for parameters
bounds = ([0, 0, 0, 0], [1, 1, 2, 1])  # Lower and upper bounds for a1, k1, a2, k2

# Fit the curve to the data
params, covariance = curve_fit(biexponential_decay, nai, nai[1], p0=(0.5, .001, 0.3, .001), bounds=bounds)  # Providing initial parameters as p0

# Extract the fitted parameters
a1_fit, k1_fit, a2_fit, k2_fit = params

# Calculate the residuals
residuals = nai[1] - biexponential_decay(nai, a1_fit, k1_fit, a2_fit, k2_fit)

# Calculate the mean squared error
mse = np.mean(residuals**2)

# Calculate the diagonal elements of the covariance matrix
cov_diag = np.diag(covariance)

# Calculate the standard error on the fitted values
std_error = np.sqrt(mse * cov_diag)

# Display the fitted parameters, R-squared value, and standard errors
print("Fitted Parameters:")
print("a1 =", a1_fit)
print("k1 =", k1_fit)
print("a2 =", a2_fit)
print("k2 =", k2_fit)
print("y0 =", 1-a1_fit-a2_fit)
print("R-squared:", 1 - mse / np.var(nai[1]))
print("Standard Errors:", std_error)

# Generate the fitted curve using the fitted parameters
y_fit = biexponential_decay(nai, a1_fit, k1_fit, a2_fit, k2_fit)

# Plot the original data and the fitted curve
plt.scatter(nai[0], nai[1], label='Original Data')
plt.plot(nai[0], y_fit, 'r-', label='Fitted Curve')
plt.legend()
plt.xlabel('nai[0]')
plt.ylabel('nai[1]')
plt.show()