引导高斯和最佳拟合数据

Question

我有以下代码，它为我的给定数据提供了最佳拟合曲线。我不知道如何引导它并使用它来找到 95% 置信区间。我的数据附在这里。

这是我的代码：

import numpy as np
import matplotlib.pyplot as plt
import pylab as plb
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])

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()

Answer 1

假设 x、y 是由高斯生成的，但存在错误：

x = np.arange(-100,100,step=0.1)

a = 2.4
x0 = 0.1
sigma = 3

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

y = gauss(x, a, x0, sigma) + 0.1 * np.random.normal(size = x.shape[0])

然后要在任何点（例如 x=0.2）生成置信区间，您可以对参数进行采样，使用该参数应用高斯函数并在该点获取 y。然后获取 10% 和 90% 分位数来绘制置信区间：

def conf90_gauss(x, popt, pcov):
    arr_y = []
    for t in range(100):
        a, mean, sigma = np.random.multivariate_normal(popt, pcov)
        y = gauss(x,a,x0,sigma)
        arr_y.append(y)

    return arr_y


conf_interval = conf90_gauss(0.2, popt, pcov)
q10, q90 = np.quantile(conf_interval, [0.1,0.9])
print("Confidence interval:", q10, q90)