Mit scipy optimieren für MLE Schätzung und Kurvenanpassung

Question

Ich generierte zufällig 1000 Datenpunkte mit den Gewichten, die ich weiß, sind wahr für die Normalverteilung. Jetzt versuche ich, die Log-Likelihood-Funktion zu minimieren, um die Werte von sig^2 und die Gewichte zu schätzen. Ich verstehe den Prozess konzeptionell, aber wenn ich versuche, ihn zu programmieren, bin ich einfach verloren.

Das ist mein Modell:

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

ich für eine Weile habe jetzt googeln und ich habe die scipy.stats.optimize.minimize Funktion ist gut für diese gelernt, aber ich kann nicht es richtig zu arbeiten. Jede Lösung, die ich versucht habe, hat für das Beispiel funktioniert, von dem ich die Lösung bekommen habe, aber ich kann es nicht auf mein Problem übertragen.

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=??)) 

habe ich gefunden, dass der Gleichungsfehler ( w) = .5 * Summe (Poly (x_n, w) - y_n)^2 relevant ist für den Fehler der Gewichtungen minimiert werden, die deshalb maximiert meine Wahrscheinlichkeit für die Gewichte, aber ich verstehe nicht, wie ich das schreiben soll ... Ich habe eine ähnliche Beziehung für sig^2 gefunden, aber habe das gleiche Problem. Kann jemand erklären, wie man das macht, um meine Kurvenanpassung zu unterstützen? Vielleicht so weit gehen, um Pseudo-Code zu posten, den ich verwenden kann?

Answer 1

Ja, die Implementierung der Likelihood-Anpassung mit minimize ist schwierig, ich verbringe viel Zeit damit. Deshalb habe ich es eingepackt. Wenn ich schamlos mein eigenes Paket symfit stecken kann, kann Ihr Problem, indem Sie so etwas wie dies gelöst werden:

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() 

Ich muss zugeben, dass ich nicht genau, Distribution verstehen, da ich nicht verstehen, die Rolle der Ihre w, aber vielleicht mit diesem Code als Beispiel, werden Sie wissen, wie Sie es anpassen.

Wenn nicht, lassen Sie mich die vollständige mathematische Gleichung Ihres Modells wissen, damit ich Ihnen weiterhelfen kann.

Für weitere Informationen überprüfen Sie die docs. (Für eine technische Beschreibung dessen, was unter der Motorhaube passiert, lesen Sie here und here.)

Answer 2

Ich denke, es gibt ein Problem mit Ihrer Einrichtung. Mit maximaler Wahrscheinlichkeit erhalten Sie die Parameter, die die Wahrscheinlichkeit der Beobachtung Ihrer Daten (bei einem bestimmten Modell) maximieren. Ihr Modell scheint zu sein:

wobei Epsilon ist N (0, Sigma).

So maximieren Sie es:

oder äquivalent Protokolle nehmen zu bekommen:

Das f in diesem Fall ist die log-Normalwahrscheinlichkeitsdichte der Du Holen Sie sich mit stats.norm.logpdf. Sie sollten dann verwenden, um einen Ausdruck zu maximieren, der die Summe von stats.norm.logpdf ist, die an jedem der i Punkte ausgewertet wird, von 1 bis zu Ihrer Stichprobengröße.

Wenn ich Sie richtig verstanden habe, fehlt Ihr Code mit einem y Vektor plus einem x Vektor! Zeigen Sie uns ein Beispiel dieser Vektoren und ich kann meine Antwort aktualisieren, um einen Beispielcode für die Schätzung von MLE mit diesem Datum aufzunehmen.