如何将curve_fit替换为scipy.optimize.minimize实现高斯拟合？

我现有一段可正常运行的代码，当前使用curve_fit方法完成高斯拟合。我希望改用scipy.optimize.minimize手动执行拟合操作，请问该方案是否可行？具体应如何实现？

原代码如下：

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

data = np.loadtxt('gaussian.dat')
x = data[:, 0]
y = data[:, 1]

n = len(x)                          
mean = sum(x*y)/n                   
sigma = sum(y*(x-mean)**2)/n        

def gauss(x,a,x0,sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))

popt,pcov = curve_fit(gauss,x,y,p0=[1,mean,sigma])

def get_conf_gaus(x: float,popt: np.ndarray, pcov: np.ndarray, n_boostrap:int = 100):

    params = np.random.multivariate_normal(popt, pcov, size = [n_boostrap])
    a = params[:,0]
    mu = params[:,1]
    sigma = params[:,2]

    bootstrap_y = gauss(0.2,a,mu,sigma)
    conf = np.quantile(bootstrap_y, [0.025,0.975])
    return conf

conf = get_conf_gaus(0.2, popt, pcov)
print(conf)

plt.plot(x,y,'b+:',label='data')
plt.plot(x,gauss(x,*popt),'ro:',label='fit')
plt.legend()
plt.title('Gaussian Fit vs Actual Data')
plt.xlabel('x-values')
plt.ylabel('y-values')
plt.show()

该方案完全可行。scipy.optimize.minimize可以通过定义损失函数（比如最小二乘损失）实现高斯拟合，核心逻辑和curve_fit一致（都是最小化残差平方和），只是需要手动构建优化目标。

实现步骤 #

• 定义损失函数：以高斯模型预测值与真实值的残差平方和作为优化目标，输入为高斯模型的参数[a, x0, sigma]。
• 初始化参数：沿用原代码基于数据均值、标准差得到的初始猜测值。
• 调用minimize：选择合适的优化器（如L-BFGS-B，适合连续参数优化）执行最小化。
• 处理置信区间：minimize不直接返回参数协方差矩阵，可通过数据bootstrap方式获取，适配原代码的置信区间计算逻辑。

完整替换后的代码 #

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

data = np.loadtxt('gaussian.dat')
x = data[:, 0]
y = data[:, 1]

n = len(x)                          
mean = sum(x*y)/n                   
sigma_init = np.sqrt(sum(y*(x-mean)**2)/n)        

def gauss(x,a,x0,sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))

# 定义最小二乘损失函数
def loss(params):
    a, x0, sigma = params
    y_pred = gauss(x, a, x0, sigma)
    return np.sum((y - y_pred)**2)

# 初始化参数
initial_guess = [1, mean, sigma_init]

# 执行最小化优化
result = minimize(loss, initial_guess, method='L-BFGS-B')
popt = result.x  # 拟合得到的最优参数

# 适配bootstrap置信区间计算（数据bootstrap方式）
def get_conf_gaus(x_val: float, popt: np.ndarray, n_bootstrap:int = 1000):
    confs = []
    for _ in range(n_bootstrap):
        # 对原始数据添加噪声后重新拟合
        y_noisy = y + np.random.normal(0, np.std(y - gauss(x, *popt)), size=len(y))
        def loss_bootstrap(params):
            a, x0, sigma = params
            return np.sum((y_noisy - gauss(x, a, x0, sigma))**2)
        res_bootstrap = minimize(loss_bootstrap, popt, method='L-BFGS-B')
        confs.append(gauss(x_val, *res_bootstrap.x))
    return np.quantile(confs, [0.025, 0.975])

conf = get_conf_gaus(0.2, popt)
print(conf)

plt.plot(x,y,'b+:',label='data')
plt.plot(x,gauss(x,*popt),'ro:',label='fit')
plt.legend()
plt.title('Gaussian Fit vs Actual Data (minimize version)')
plt.xlabel('x-values')
plt.ylabel('y-values')
plt.show()

关键说明 #

• 损失函数：loss函数计算残差平方和，和curve_fit默认优化目标一致，因此拟合结果和原代码几乎相同。
• 优化器选择：L-BFGS-B适合连续参数优化，若无需梯度信息，也可选择Nelder-Mead等无梯度优化器。
• 置信区间：改用数据bootstrap方式（对原始数据添加噪声后重新拟合），比原代码的参数bootstrap更严谨，避免依赖协方差矩阵的假设。

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