正态分布的MLE拟合

Question

我试图采用this thread中提出的这个解来确定一个简单正态分布的参数。尽管修改是很小的(基于维基百科)，但结果却很差。有什么不对的地方吗？

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:数组(2.27373675e-05，2.27373675e-05)

信息：“会聚: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样本是您的实际目标，则可以使用曲线拟合函数(如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，则必须相应地更改绘图代码。)