对数正态分布的可能性

Question

我试图使用最小化函数形式scipy找到正态分布的mu和sigma的最大似然估计。然而，最小化返回平均值的预期值，但西格玛的估计远非真实的西格玛。

我定义函数llnorm返回正态分布的负对数似然，然后从正态分布创建随机样本，平均值为150，标准差为10，然后使用optimize我试图找到MLE。

import numpy as np
import math
import scipy.optimize as optimize

def llnorm(par, data):
    n = len(data)
    mu, sigma = par
    ll = -np.sum(-n/2 * math.log(2*math.pi*(sigma**2)) - ((data-mu)**2)/(2 * (sigma**2)))
    return ll

data = 10 * np.random.randn(100) + 150

result = optimize.minimize(llnorm, [150,10], args = (data))

即使数据的均值接近150且std接近10，最优化也会使估计的sigma值（接近0）返回得更小。

Answer 1

你的数学略有偏差：

ll = n*math.log(2*math.pi*(sigma**2))/2 + np.sum(((data-mu)**2)/(2 * (sigma**2)))

要么

ll = np.sum(math.log(2*math.pi*(sigma**2))/2 + ((data-mu)**2)/(2 * (sigma**2)))

首先我取消了-（不是问题），但最重要的是你要保持总和中的常数项并且不要将它乘以n，或者你把它取出并乘以n，...但不是两者都在同一时间。

Answer 2

np.random.randn创建随机高斯分布，方差为1（docs here）。由于你的目标是使用std为10的分布，你需要用10 * 10相乘

import numpy as np
import math
import scipy.optimize as optimize

def llnorm(par, data):
    n = len(data)
    mu, sigma = par
    ll = -np.sum(-n/2 * math.log(2*math.pi*(sigma**2)) - ((data-mu)**2)/(2 * (sigma**2)))
    return ll

data = 10 * 10 * np.random.randn(100) + 150 

result = optimize.minimize(llnorm, [150,10], args = (data))
print(result)

这给了我：

      fun: 36328.17002555693
 hess_inv: array([[ 0.96235834, -0.32116447],
       [-0.32116447,  0.10879383]])
      jac: array([0., 0.])
  message: 'Optimization terminated successfully.'
     nfev: 44
      nit: 8
     njev: 11
   status: 0
  success: True
        x: array([166.27014352,   9.15113937])

编辑：似乎~9的输出纯属巧合。还有其他东西需要调查