Wie kann man MLE-Schätzer in Python numerisch schätzen, wenn die Gradienten sehr klein sind und weit von der optimalen Lösung entfernt sind?

Ich untersuche, wie man einen Datensatz unter Verwendung von Normalverteilungen modelliert, wobei sowohl Mittelwert als auch Varianz als lineare Funktionen unabhängiger Variablen definiert sind.

So etwas wie N ~ (f (x), g (x)).

Ich generiere eine Zufallsstichprobe wie diese:

def draw(x):
    return norm(5 * x + 2, 3 *x + 4).rvs(1)[0]

Ich möchte also 5, 2 und 4 als Parameter für meine Verteilung abrufen.

Ich generiere meine Probe:

smp = np.zeros ((100,2))

for i in range(0, len(smp)):
    smp[i][0] = i
    smp[i][1] = draw(i)

Die Wahrscheinlichkeitsfunktion ist:

def lh(p):
    p_loc_b0 = p[0]
    p_loc_b1 = p[1]
    p_scl_b0 = p[2]
    p_scl_b1 = p[3]

    l = 1
    for i in range(0, len(smp)):
        x = smp[i][0]
        y = smp[i][1]
        l = l * norm(p_loc_b0 + p_loc_b1 * x, p_scl_b0 + p_scl_b1 * x).pdf(y)

    return -l

Die Parameter für die im Modell verwendeten linearen Funktionen sind also im p 4 -variablen Vektor angegeben.

Mit scipy.optimize kann ich die MLE-Parameter mit einem extrem niedrigen xtol lösen und die Lösung bereits als Ausgangspunkt angeben:

fmin(lh, x0=[2,5,3,4], xtol=1e-35)

Was nicht gut funktioniert:

Warning: Maximum number of function evaluations has been exceeded.
array([ 3.27491346,  4.69237042,  5.70317719,  3.30395462])

Das Erhöhen des xtol auf höhere Werte bringt nichts.

Also versuche ich, eine Startlösung zu verwenden, die weit von der tatsächlichen Lösung entfernt ist:

>>> fmin(lh, x0=[1,1,1,1], xtol=1e-8)
Optimization terminated successfully.
         Current function value: -0.000000
         Iterations: 24
         Function evaluations: 143
array([ 1.,  1.,  1.,  1.])

Was mich zum Nachdenken bringt:

PDF-Dateien sind größtenteils um den Mittelwert gruppiert und weisen nur wenige Standardabweichungen vom Mittelwert auf, die für numerische Methoden nicht zu gut sein dürfen.

Wie führt man diese Art der numerischen Schätzung in Funktionen durch, bei denen der Gradient sehr nahe bei Null von der Lösung entfernt ist?

Es gibt mehrere Gründe, warum Sie fehlerhafte Ergebnisse erhalten. Zunächst sollten Sie die Verwendung der Protokollwahrscheinlichkeit anstelle der Wahrscheinlichkeit in Betracht ziehen. Es gibt numerische Probleme beim Multiplizieren vieler kleiner Zahlen (stellen Sie sich vor, wenn Sie Millionen von Stichproben hätten, müssten Sie Millionen kleiner Zahlen für die lhd multiplizieren). Auch das Aufnehmen von Verläufen für Optimierungsmethoden, die Verläufe erfordern, ist häufig einfacher, wenn Sie sich mit der Protokollwahrscheinlichkeit befassen. Im Allgemeinen ist es gut, ein Ziel zu haben, das bei Optimierungsproblemen eher eine Summe als ein Produkt von Variablen ist.

Zweitens verwendet fmin den Nelder-Mead-Simplex-Algorithmus, für den laut Scipy-Dokumentation keine Konvergenzgarantien bestehen . Dies bedeutet, dass die Konvergenz völlig zufällig ist und Sie nicht erwarten sollten, Parameter in der Nähe der Originale zu finden. Um dies zu umgehen, würde ich Ihnen empfehlen, eine gradientenbasierte Methode wie stochastischen Gradientenabstieg oder BFGS zu verwenden. Da Sie das generative Modell kennen (rvs sind Gauß-verteilt), können Sie die Wahrscheinlichkeit und die Protokollwahrscheinlichkeit wie folgt schreiben:

Wobei a, b, c und d Ihre Modellparameter 5,2,3 bzw. 4 sind. Nehmen Sie dann den Gradienten in Bezug auf [a, b, c, d] und geben Sie ihn in die Primäreingabe von fmin_bfgs ein. Beachten Sie, dass aufgrund der unterschiedlichen Varianz das, was nur durch lineare Regression gelöst werden könnte, jetzt ein schlimmeres Problem ist.

Schließlich möchten Sie möglicherweise auch die verallgemeinerten kleinsten Quadrate auf http://en.wikipedia.org/wiki/Linear_regression#Least-squares_estimation_and_related_techniques und http://en.wikipedia.org/wiki/Heteroscedasticity überprüfen , die über Ihr Problem und Ihre Probleme sprechen bieten mehrere verfügbare Lösungen.

Viel Glück!