Блог

Я случайным образом сгенерировал 1000 точек данных, используя веса, которые, как я знаю, верны для нормального распределения. Теперь я пытаюсь минимизировать функцию логарифмического правдоподобия для оценки значений sig ^ 2 и весов. Я вроде бы понимаю процесс концептуально, но когда я пытаюсь его закодировать, я просто теряюсь.

Это моя модель:

 p(y|x, w, sig^2) = N(y|w0 w1x ... wnx^n, sig^2)
  

Я некоторое время гуглил и узнал, что функция scipy.stats.optimize.minimize хороша для этого, но я не могу заставить ее работать правильно. Каждое решение, которое я пробовал, работало для примера, из которого я получил решение, но я не могу экстраполировать его на свою проблему.

 x = np.linspace(0, 1000, num=1000)
data = []
for y in x:
        data.append(np.polyval([.5, 1, 3], y))

#plot to confirm I do have a normal distribution...
data.sort()
pdf = stats.norm.pdf(data, np.mean(data), np.std(data))
plt.plot(test, pdf)
plt.show()

#This is where I am stuck.
logLik = -np.sum(stats.norm.logpdf(data, loc=??, scale=??))
  

Я обнаружил, что ошибка уравнения (w) = .5 * sum(poly(x_n, w) — y_n) ^ 2 важна для минимизации ошибки весов, что, следовательно, максимизирует мою вероятность для весов, но я не понимаю, как это закодировать… Я нашел аналогичную взаимосвязь для sig ^ 2, но столкнулся с той же проблемой. Кто-нибудь может пояснить, как это сделать, чтобы помочь моей подгонке кривой? Может быть, зайти так далеко, чтобы опубликовать псевдокод, который я могу использовать?

Да, реализовать подгонку правдоподобия с помощью minimize сложно, я трачу на это много времени. Вот почему я это завернул. Если я могу без зазрения совести подключить свой собственный пакет symfit , вашу проблему можно решить, выполнив что-то вроде этого:

 from symfit import Parameter, Variable, Likelihood, exp
import numpy as np

# Define the model for an exponential distribution
beta = Parameter()
x = Variable()
model = (1 / beta) * exp(-x / beta)

# Draw 100 samples from an exponential distribution with beta=5.5
data = np.random.exponential(5.5, 100)

# Do the fitting!
fit = Likelihood(model, data)
fit_result = fit.execute()
  

Должен признать, я не совсем понимаю ваш дистрибутив, поскольку не понимаю вашей роли w , но, возможно, на примере этого кода вы поймете, как его адаптировать.

Если нет, сообщите мне полное математическое уравнение вашей модели, чтобы я мог помочь вам в дальнейшем.

Для получения дополнительной информации ознакомьтесь с документами. (Более техническое описание того, что происходит под капотом, читайте здесь и здесь.)

Я думаю, что есть проблема с вашей настройкой. С максимальной вероятностью вы получаете параметры, которые максимизируют вероятность наблюдения ваших данных (при определенной модели). Ваша модель, похоже,:

где epsilon равно N (0, сигма).

Таким образом, вы максимизируете его:

или, что эквивалентно, возьмите журналы, чтобы получить:

F в данном случае является логарифмически нормальной функцией плотности вероятности, которую вы можете получить с помощью stats.norm.logpdf . Затем вы должны использовать scipy.minimize для максимизации выражения, которое будет представлять собой суммирование stats.norm.logpdf , вычисленных в каждой из точек i, от 1 до вашего размера выборки.

Если я вас правильно понял, в вашем коде отсутствует y вектор плюс x вектор! Покажите нам образец этих векторов, и я смогу обновить свой ответ, включив пример кода для оценки MLE с этой датой.