Лог вероятности нормального распределения - PullRequest
Новая
       10

Лог вероятности нормального распределения

0 голосов
/ 04 апреля 2019

Я пытаюсь найти оценку максимального правдоподобия мю и сигмы из нормального распределения, используя функцию минимизации scipy. Однако при минимизации возвращается ожидаемое значение среднего значения, но оценка сигмы далека от реальной сигмы.

Я определяю функцию llnorm, которая возвращает отрицательное логарифмическое правдоподобие нормального распределения, затем создаю случайную выборку из нормального распределения со средним значением 150 и стандартным отклонением 10, затем с помощью оптимизации пытаюсь найти 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 и стандартное значение близко к 10, оптимизация возвращает намного меньшее значение оценочной сигмы (близко к 0).

Ответы [ 2 ]

1 голос
/ 04 апреля 2019

Ваша математика немного не в порядке: 

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, ... но не оба одновременно.

0 голосов
/ 04 апреля 2019

np.random.randn создает случайное распределение Гаусса с дисперсией 1 ( документы здесь ).Поскольку вы стремитесь получить распределение со стандартным значением 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 является чисто случайным.Что-то еще нужно исследовать

...