Как мне ограничить сумму нескольких параметров подгонки до 1 при подгонке биэкспоненты в Python? ⇐ Python

[Anonymous]

У меня есть некоторые данные о затухании, которые выглядят так, будто они могут быть биэкспоненциальными, поэтому я хочу подогнать их к уравнению y(t) = y0 + a1e^(-k1t) + a2< em>e^(-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()


Подробнее здесь: https://stackoverflow.com/questions/785 ... itting-a-b