Использование scipy оптимизация для оценки MLE и подгонки кривой

Question

Я произвольно сгенерировал 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=??)) 

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

Answer 1

Да, реализация правдоподобия с 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, но, возможно, с этим кодом в качестве примера вы узнаете, как его адаптировать.

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

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

Answer 2

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

где эпсилон является N (0, сигма).

Так вы увеличиваете его:

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

Диафрагменной в этом случае есть вероятность функция плотности Логнормальной которой вы можете получить с stats.norm.logpdf. Затем вы должны использовать , чтобы максимизировать выражение, которое будет суммировать stats.norm.logpdf, оцененное в каждой из точек i, от 1 до вашего размера выборки.

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