如何将多高斯拟合函数funcn传入scipy curve_fit并设置bounds

问题背景 #

已知X为多个高斯函数的和，现有双高斯拟合代码（当Cluster=2时可正常运行）：

双高斯函数定义：

def func2(x, a1, a2, m1, m2, s1, s2):
    return a1*np.exp(-((x-m1)/s1)**2) + a2*np.exp(-((x-m2)/s2)**2)

配合curve_fit的调用代码：

# 设置参数范围
AmpMin = 0.01
AmpMax = 1000
CenMin = min(x)
CenMax = max(x)
WidMin = 0.1
WidMax = 100

popt, pcov = curve_fit(func2, x, y, 
                       bounds=([AmpMin,AmpMin, CenMin,CenMin, WidMin,WidMin],
                               [AmpMax,AmpMax, CenMax,CenMax, WidMax,WidMax]))

现在需要将上述代码泛化为支持指定Cluster数量的多高斯拟合，但原有的funcn写法无法直接传入curve_fit：

def fun1(x,a,m,s):
    return a*np.exp(-((x-m)/s)**2)

def funcn(x, as[],ms[],ss[]):
    accum = fun1(x,as[0],ms[0],ss[0])
    for i in range(1,len(as)):
        accum = accum + fun1(x,as[i],ms[i],ss[i])
    return accum

需解决：如何调整funcn使其适配scipy.optimize.curve_fit，并正确设置bounds参数？

解决方案 #

核心问题说明 #

scipy.optimize.curve_fit要求拟合函数的所有待优化参数必须是扁平的位置参数，不能直接传入列表类型的参数，因此需要重构多高斯函数的参数结构。

1. 重构多高斯拟合函数 #

将所有高斯分量的振幅(a)、中心(m)、宽度(s)作为扁平的位置参数传入，在函数内部按每组3个参数（对应一个高斯分量）拆分处理：

import numpy as np
from scipy.optimize import curve_fit

def single_gaussian(x, a, m, s):
    """单个高斯函数"""
    return a * np.exp(-((x - m) / s) ** 2)

def multi_gaussian(x, *params):
    """支持任意数量高斯分量的拟合函数
    params参数格式：(a1, m1, s1, a2, m2, s2, ..., an, mn, sn)
    """
    n_clusters = len(params) // 3
    if len(params) % 3 != 0:
        raise ValueError("参数数量必须是3的倍数（每个高斯分量对应a、m、s三个参数）")
    
    y = np.zeros_like(x)
    for i in range(n_clusters):
        a = params[3 * i]
        m = params[3 * i + 1]
        s = params[3 * i + 2]
        y += single_gaussian(x, a, m, s)
    return y

2. 生成适配的bounds参数 #

根据指定的Cluster数量，按(a_min, m_min, s_min)重复n次生成下界，(a_max, m_max, s_max)重复n次生成上界：

# 定义参数范围
AmpMin = 0.01
AmpMax = 1000
CenMin = x.min()
CenMax = x.max()
WidMin = 0.1
WidMax = 100

# 指定高斯分量数量
n_clusters = 3  # 可根据需求修改为任意正整数

# 生成bounds
lower_bounds = [AmpMin, CenMin, WidMin] * n_clusters
upper_bounds = [AmpMax, CenMax, WidMax] * n_clusters

3. 调用curve_fit（需注意初始值p0） #

多峰拟合的初始猜测值p0非常关键，合理的初始值能避免拟合陷入局部最优。以下是示例初始值生成逻辑（可根据实际数据调整）：

# 生成初始猜测值p0
p0 = []
for i in range(n_clusters):
    # 中心位置在数据范围内均匀分布
    cen = CenMin + (CenMax - CenMin) * (i + 0.5) / n_clusters
    # 振幅初始为数据最大值的1/n
    amp = y.max() / n_clusters
    # 宽度初始为数据范围的1/(2n)
    wid = (CenMax - CenMin) / (2 * n_clusters)
    p0.extend([amp, cen, wid])

# 执行拟合
popt, pcov = curve_fit(multi_gaussian, x, y, p0=p0, bounds=(lower_bounds, upper_bounds))

# 解析拟合结果：拆分每个高斯分量的参数
fit_results = []
for i in range(n_clusters):
    a = popt[3 * i]
    m = popt[3 * i + 1]
    s = popt[3 * i + 2]
    fit_results.append({"振幅": a, "中心": m, "宽度": s})

可选参数顺序调整 #

如果偏好先传入所有振幅、再所有中心、再所有宽度的参数顺序，只需修改函数内部的参数拆分逻辑，并同步调整bounds和p0的生成顺序即可：

def multi_gaussian(x, *params):
    n_clusters = len(params) // 3
    a_list = params[:n_clusters]
    m_list = params[n_clusters:2*n_clusters]
    s_list = params[2*n_clusters:]
    
    y = np.zeros_like(x)
    for a, m, s in zip(a_list, m_list, s_list):
        y += single_gaussian(x, a, m, s)
    return y

# 对应调整bounds
lower_bounds = [AmpMin]*n_clusters + [CenMin]*n_clusters + [WidMin]*n_clusters
upper_bounds = [AmpMax]*n_clusters + [CenMax]*n_clusters + [WidMax]*n_clusters

# p0也同步调整
p0 = [y.max()/n_clusters]*n_clusters + \
     [CenMin + (CenMax-CenMin)*(i+0.5)/n_clusters for i in range(n_clusters)] + \
     [(CenMax-CenMin)/(2*n_clusters)]*n_clusters

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