我试图采用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:array([2.27373675e-05,2.27373675e-05])
消息:'CONVERGENCE:REL_REDUCTION_OF_F _ <= _ FACTR * EPSMCH'
成功:是的
x:数组([10.45000245,5.48475283])
最大似然法用于将分布的参数拟合到一组值,这些值据称是来自该分布的随机样本。在你的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
作为样本,则必须相应地更改绘图代码。)