高斯拟合曲线的Bootstrapping实现及95%置信区间求解问询

方法概述 #

Bootstrapping的核心思路是通过有放回地重复采样原始数据，每次采样后重新拟合模型，收集所有拟合得到的参数，最终通过参数分布的分位数来确定置信区间。针对你的高斯拟合任务，我们将：

1. 多次重采样(x,y)数据对
2. 每次重采样后重新拟合高斯曲线
3. 统计所有拟合参数的分布，计算95%置信区间（取2.5%和97.5%分位数）
4. 绘制拟合曲线的置信带以直观展示不确定性

修改后的完整代码 #

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])
print(f"原始拟合参数: a={popt[0]:.4f}, x0={popt[1]:.4f}, sigma={popt[2]:.4f}")

# -------------------------- Bootstrapping 核心部分 --------------------------
# 设置bootstrap迭代次数（建议至少1000次，次数越多结果越稳定）
n_bootstrap = 1000

# 存储每次拟合的参数（a, x0, sigma）
boot_params = np.zeros((n_bootstrap, 3))

for i in range(n_bootstrap):
    # 有放回地采样数据索引
    idx = np.random.choice(n, size=n, replace=True)
    x_boot = x[idx]
    y_boot = y[idx]
    
    # 对重采样数据进行拟合（使用原始拟合参数作为初始猜测提升收敛性）
    try:
        popt_boot, _ = curve_fit(gauss, x_boot, y_boot, p0=popt)
        boot_params[i] = popt_boot
    except RuntimeError:
        # 若某次拟合失败（如重采样后数据分布极端），跳过该次迭代
        print(f"第{i}次拟合失败，跳过")
        continue

# 过滤掉未成功拟合的无效参数
boot_params = boot_params[~np.all(boot_params == 0, axis=1)]

# 计算95%置信区间
confidence_interval = np.percentile(boot_params, [2.5, 97.5], axis=0)
print("\n95%置信区间:")
print(f"a: [{confidence_interval[0,0]:.4f}, {confidence_interval[1,0]:.4f}]")
print(f"x0: [{confidence_interval[0,1]:.4f}, {confidence_interval[1,1]:.4f}]")
print(f"sigma: [{confidence_interval[0,2]:.4f}, {confidence_interval[1,2]:.4f}]")

# -------------------------- 可视化部分 --------------------------
plt.figure(figsize=(10,6))

# 绘制原始数据和原始拟合曲线
plt.plot(x, y, 'b+:', label='原始数据')
plt.plot(x, gauss(x, *popt), 'r-', label='原始拟合曲线', linewidth=2)

# 绘制部分bootstrap拟合曲线（展示拟合结果的变异性）
for params in boot_params[:50]:  # 仅绘制50条避免图像杂乱
    plt.plot(x, gauss(x, *params), 'gray', alpha=0.1)

# 计算并绘制拟合曲线的置信带
x_grid = np.linspace(min(x), max(x), 100)
boot_curves = np.array([gauss(x_grid, *params) for params in boot_params])
lower_band = np.percentile(boot_curves, 2.5, axis=0)
upper_band = np.percentile(boot_curves, 97.5, axis=0)
plt.fill_between(x_grid, lower_band, upper_band, color='pink', alpha=0.3, label='95%置信带')

plt.legend()
plt.title('高斯拟合与Bootstrap置信区间')
plt.xlabel('x-values')
plt.ylabel('y-values')
plt.show()

关键细节说明 #

• 重采样逻辑：采用有放回采样原始(x,y)数据对，确保每次重采样的数据集与原始数据规模一致，符合bootstrapping的核心要求。
• 拟合稳定性：加入try-except块处理拟合失败的异常情况，避免因极端重采样数据导致程序中断。
• 置信区间计算：通过np.percentile直接获取参数分布的2.5%和97.5%分位数，对应95%置信区间，无需假设参数服从特定分布。
• 可视化增强：灰色细线条展示部分bootstrap拟合结果的变异性，粉色填充区域为拟合曲线的95%置信带，直观呈现拟合的不确定性范围。

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