有没有办法从 scipy.stats.norm.fit 中获取拟合参数的错误？

Question

我有一些数据，我使用 scipy.stats.normal 对象拟合函数拟合正态分布，如下所示：

import numpy as np                                                                                                                                                                                                                       
import matplotlib.pyplot as plt                                                                                                                                                                                                          
from scipy.stats import norm                                                                                                                                                                                                             
import matplotlib.mlab as mlab                                                                                                                                                                                                           

x = np.random.normal(size=50000)                                                                                                                                                                                                         

fig, ax = plt.subplots()                                                                                                                                                                                                                 

nbins = 75                                                                                                                                                                                                                               
mu, sigma = norm.fit(x)                                                                                                                                                                                                                  
n, bins, patches = ax.hist(x,nbins,normed=1,facecolor = 'grey', alpha = 0.5, label='before');                                                                                                                                            
y0 = mlab.normpdf(bins, mu, sigma) # Line of best fit                                                                                                                                                                                    
ax.plot(bins,y0,'k--',linewidth = 2, label='fit before')                                                                                                                                                                                 
ax.set_title('$\mu$={}, $\sigma$={}'.format(mu, sigma))                                                                                                                                                                                  

plt.show()                                                                                                                                                                                                                               

我现在想提取拟合的 mu 和 sigma 值中的不确定性/误差。我该怎么办？

Answer 1

您可以使用scipy.optimize.curve_fit： 该方法不仅返回估计的最优值 参数，还有对应的协方差矩阵：

popt : 数组

参数的最优值使得残差平方和 f(xdata, *popt) - ydata 被最小化

pcov：二维数组

popt 的估计协方差。对角线提供参数估计的方差。要计算参数的一个标准差误差，请使用 perr = np.sqrt(np.diag(pcov))。

如上所述，sigma 参数如何影响估计的协方差取决于 absolute_sigma 参数。

如果解决方案中的雅可比矩阵没有满秩，则 'lm' 方法返回一个填充有 np.inf 的矩阵，另一方面，'trf' 和 'dogbox' 方法使用 Moore-Penrose 伪逆计算协方差矩阵。

您可以从协方差矩阵的对角元素的平方根计算参数的标准偏差误差如下：

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

x = np.random.normal(size=50000)
fig, ax = plt.subplots() 
nbins = 75
n, bins, patches = ax.hist(x,nbins, density=True, facecolor = 'grey', alpha = 0.5, label='before'); 

centers = (0.5*(bins[1:]+bins[:-1]))
pars, cov = curve_fit(lambda x, mu, sig : norm.pdf(x, loc=mu, scale=sig), centers, n, p0=[0,1])

ax.plot(centers, norm.pdf(centers,*pars), 'k--',linewidth = 2, label='fit before') 
ax.set_title('$\mu={:.4f}\pm{:.4f}$, $\sigma={:.4f}\pm{:.4f}$'.format(pars[0],np.sqrt(cov[0,0]), pars[1], np.sqrt(cov[1,1 ])))

plt.show()

这导致以下情节：

Answer 2

另请参阅 lmfit (https://github.com/lmfit/lmfit-py)，它提供了更简单的界面并报告拟合变量的不确定性。要将数据拟合到正态分布，请参阅http://lmfit.github.io/lmfit-py/builtin_models.html#example-1-fit-peak-data-to-gaussian-lorentzian-and-voigt-profiles

并使用类似的东西

from lmfit.models import GaussianModel

model = GaussianModel()

# create parameters with initial guesses:
params = model.make_params(center=9, amplitude=40, sigma=1)  

result = model.fit(ydata, params, x=xdata)
print(result.fit_report())

报告将包括 1-sigma 错误，例如

[[Variables]]
    sigma:       1.23218358 +/- 0.007374 (0.60%) (init= 1.0)
    center:      9.24277047 +/- 0.007374 (0.08%) (init= 9.0)
    amplitude:   30.3135620 +/- 0.157126 (0.52%) (init= 40.0)
    fwhm:        2.90157055 +/- 0.017366 (0.60%)  == '2.3548200*sigma'
    height:      9.81457817 +/- 0.050872 (0.52%)  == '0.3989423*amplitude/max(1.e-15, sigma)'