如何在Python中为光变曲线拟合高斯上升+指数衰减模型？

我需要给光变曲线数据拟合一种峰值前高斯上升、峰值后指数衰减的模型，目标曲线为蓝色观测曲线。现有拟合高斯模型的初始代码，想实现：限制高斯函数仅作用于峰值前数据，指数/幂律函数从峰值开始拟合后续数据。

初始代码 #

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
 
# 指数衰减函数
def exponential(x, a, b):
    return a*np.exp(-b*x)
 
# 幂律衰减函数
def power_law(x, a, b):
    return a*np.power(-x, b)

# 高斯上升函数
def gaussian(x, a, b, c):
    return a*np.exp(-np.power(x - b, 2)/(2*np.power(c, 2)))

# 数据准备（假设t_slice是已加载的DataFrame）
xData = t_slice['time']
yData = t_slice['flux']/10**38 # 归一化处理

# 绘图
plt.plot(xData, yData, 'b.', label='SED')
 
# 拟合高斯模型
pars, cov = curve_fit(gaussian, xdata=xData, ydata=yData, p0=[0, 0, 0], bounds=(-np.inf, np.inf))
stdevs = np.sqrt(np.diag(cov))

print("高斯参数:",pars)
print("参数标准差:",stdevs)
 
plt.scatter(xData, gaussian(xData, *pars), linestyle='--', linewidth=2, color='black')

plt.ylim(0, 1000)
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

解决方案：构建分段复合拟合模型 #

核心思路是创建一个分段函数，根据输入的时间点判断处于峰值前还是后，分别调用高斯或指数/幂律函数，同时保证拟合的稳定性和连续性。

方案1：手动确定峰值位置（适合峰值明显的情况） #

先从观测数据中提取峰值对应的时间点，再基于该点拆分数据进行分段拟合：

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

# 1. 提取数据中的峰值位置
peak_idx = np.argmax(yData)
peak_x = xData.iloc[peak_idx] if isinstance(xData, pd.Series) else xData[peak_idx]

# 2. 定义分段复合函数：峰值前高斯，峰值后指数衰减
def piecewise_model(x, a_gauss, b_gauss, c_gauss, a_exp, b_exp):
    mask = x <= peak_x
    y = np.zeros_like(x)
    # 峰值前：高斯上升
    y[mask] = a_gauss * np.exp(-np.power(x[mask] - b_gauss, 2)/(2*np.power(c_gauss, 2)))
    # 峰值后：从峰值点开始的指数衰减（用x-peak_x保证衰减起始于峰值）
    y[~mask] = a_exp * np.exp(-b_exp * (x[~mask] - peak_x))
    return y

# 3. 设置初始参数：高斯部分复用之前的拟合结果，指数部分初始值设为峰值处的流量
p0 = [pars[0], pars[1], pars[2], yData[peak_idx], 0.001]
# 4. 设置参数边界：避免无意义值（如标准差为负、衰减系数为负）
bounds = ([-np.inf, -np.inf, 1e-3, 0, 0], [np.inf, np.inf, np.inf, np.inf, np.inf])
pars_piecewise, cov_piecewise = curve_fit(piecewise_model, xData, yData, p0=p0, bounds=bounds)

# 5. 生成拟合曲线并绘图
xFit = np.linspace(xData.min(), xData.max(), 1000)
yFit = piecewise_model(xFit, *pars_piecewise)

plt.plot(xData, yData, 'b.', label='观测数据')
plt.plot(xFit, yFit, 'r--', linewidth=2, label='分段拟合曲线')
plt.axvline(x=peak_x, color='gray', linestyle=':', label='峰值位置')
plt.ylim(0, 1000)
plt.xlabel('时间')
plt.ylabel('归一化流量')
plt.legend()
plt.show()

方案2：将峰值位置作为拟合参数（更灵活） #

无需手动找峰值，让优化器自动求解最优峰值位置及两段函数的参数：

def piecewise_model_with_peak(x, peak_x, a_gauss, b_gauss, c_gauss, a_exp, b_exp):
    mask = x <= peak_x
    y = np.zeros_like(x)
    # 峰值前高斯上升
    y[mask] = a_gauss * np.exp(-np.power(x[mask] - b_gauss, 2)/(2*np.power(c_gauss, 2)))
    # 峰值后指数衰减
    y[~mask] = a_exp * np.exp(-b_exp * (x[~mask] - peak_x))
    return y

# 初始参数：峰值位置设为数据中流量最大的时间点
initial_peak_x = xData.iloc[np.argmax(yData)] if isinstance(xData, pd.Series) else xData[np.argmax(yData)]
p0 = [initial_peak_x, pars[0], pars[1], pars[2], yData[np.argmax(yData)], 0.001]
# 参数边界：峰值位置必须在数据的时间范围内
bounds = ([xData.min(), -np.inf, -np.inf, 1e-3, 0, 0], [xData.max(), np.inf, np.inf, np.inf, np.inf, np.inf])
pars_piecewise, cov_piecewise = curve_fit(piecewise_model_with_peak, xData, yData, p0=p0, bounds=bounds)

# 绘图
xFit = np.linspace(xData.min(), xData.max(), 1000)
yFit = piecewise_model_with_peak(xFit, *pars_piecewise)

plt.plot(xData, yData, 'b.', label='观测数据')
plt.plot(xFit, yFit, 'r--', linewidth=2, label='分段拟合曲线')
plt.axvline(x=pars_piecewise[0], color='gray', linestyle=':', label='拟合峰值位置')
plt.ylim(0, 1000)
plt.xlabel('时间')
plt.ylabel('归一化流量')
plt.legend()
plt.show()

替换为幂律衰减 #

若需要用幂律替代指数衰减，仅需修改分段函数的峰值后部分：

# 定义幂律衰减函数（x-peak_x为正数，无需负号）
def power_law_decay(x, a_pl, b_pl):
    return a_pl * np.power(x, b_pl)

# 修改分段模型的峰值后逻辑
y[~mask] = a_pl * np.power(x[~mask] - peak_x, b_pl)
# 注意调整初始参数和边界：b_pl通常设为负数（保证衰减趋势）

关键注意事项 #

• 初始参数设置：尽量让初始值接近真实情况（比如高斯参数复用之前的拟合结果），否则拟合可能不收敛
• 参数边界：合理设置边界避免无意义解（如高斯标准差不能为0、衰减系数必须为正）
• 连续性保证：若需要峰值点处导数连续，可添加约束条件，但会增加拟合复杂度，一般函数值连续即可满足需求

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