Scipy 正态分布的 MLE 拟合

Question

我试图采用中提出的这个解决方案来确定简单正态分布的参数。即使修改很小（基于维基百科），结果也很差。有什么地方出错的建议吗？

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

def gaussian(x, mu, sig):
    return 1./(math.sqrt(2.*math.pi)*sig)*np.exp(-np.power((x - mu)/sig, 2.)/2)

def lik(parameters):

    mu    = parameters[0]
    sigma = parameters[1]    
    n     = len(x)  
    L     = n/2.0 * np.log(2 * np.pi) + n/2.0 * math.log(sigma **2 ) + 1/(2*sigma**2) * sum([(x_ - mu)**2 for x_ in x ])

    return L



mu0    = 10
sigma0 = 2

x = np.arange(1,20, 0.1)
y = gaussian(x, mu0, sigma0)


lik_model = minimize(lik, np.array([5,5]), method='L-BFGS-B')


mu    = lik_model['x'][0]
sigma = lik_model['x'][1]

print lik_model

plt.plot(x, gaussian(x, mu, sigma), label = 'fit')
plt.plot(x, y, label = 'data')
plt.legend()

拟合输出：

jac: array([2.27373675e-05, 2.27373675e-05])

留言：'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'

成功：正确

x: 数组([10.45000245, 5.48475283])

Answer 1

最大似然法用于将分布参数拟合到一组值，据称这些值是来自该分布的随机样本。在您的 lik 函数中，您使用 x 来保存样本，但是 x 是您已设置为 x = np.arange(1,20, 0.1) 的全局变量。那绝对不是来自正态分布的随机样本。

因为您使用的是正态分布，所以您可以使用最大似然估计的已知公式来检查您的计算：mu 是样本均值，sigma 是样本标准差：

In [17]: x.mean() Out[17]: 10.450000000000006 In [18]: x.std() Out[18]: 5.484751589634671

这些值与您调用 minimize 的结果非常接近，因此看起来您的代码正在运行。

要修改您的代码以按照您预期的方式使用 MLE，x 应该是据称是来自正态分布的随机样本的值的集合。请注意，您的数组 y 不是这样的样本。它是网格上概率密度函数 (PDF) 的值。如果将分布拟合到 PDF 样本是您的实际目标，您可以使用 curve-fitting 函数，例如 scipy.optimize.curve_fit。如果将正态分布参数拟合到随机样本实际上是您想要做的，那么为了测试您的代码，您应该使用来自具有已知参数的分布的相当大的样本作为输入。在这种情况下，你可以做

x = np.random.normal(loc=mu0, scale=sigma0, size=20)

当我在你的代码中使用这样的 x 时，我得到

In [20]: lik_model.x Out[20]: array([ 9.5760996 , 2.01946582])

正如预期的那样，解决方案中的值大约为 10 和 2。

（如果您像我一样使用 x 作为样本，则必须更改您的相应地绘制代码。）

Scipy 正态分布的 MLE 拟合

Scipy MLE fit of a normal distribution

python

gaussian

curve-fitting

scipy