eine MLE für ein Paar von Verteilungen in Python

Question

Ok Producing, so hat mein aktueller Kurvenanpassungscode einen Schritt, den scipy.stats verwendet die richtige Verteilung auf den Daten, um zu bestimmen, basierend

distributions = [st.laplace, st.norm, st.expon, st.dweibull, st.invweibull, st.lognorm, st.uniform] 
mles = [] 

for distribution in distributions: 
    pars = distribution.fit(data) 
    mle = distribution.nnlf(pars, data) 
    mles.append(mle) 

results = [(distribution.name, mle) for distribution, mle in zip(distributions, mles)] 

for dist in sorted(zip(distributions, mles), key=lambda d: d[1]): 
    print dist 
best_fit = sorted(zip(distributions, mles), key=lambda d: d[1])[0] 
print 'Best fit reached using {}, MLE value: {}'.format(best_fit[0].name, best_fit[1])   


print [mod[0].name for mod in sorted(zip(distributions, mles), key=lambda d: d[1])] 

wo Daten ein Liste der numerischen Werte. Dies funktioniert bisher sehr gut für die Anpassung unimodaler Verteilungen, bestätigt in einem Skript, das zufällig Werte aus zufälligen Verteilungen erzeugt und curve_fit verwendet, um die Parameter neu zu bestimmen.

Nun würde Ich mag den Code in der Lage machen bimodale Verteilungen zu handhaben, wie im folgenden Beispiel:

Ist es möglich, eine MLE für ein Paar Modelle zu erhalten aus scipy.stats, um festzustellen, ob ein bestimmtes Verteilungspaar für die Daten geeignet ist ?, so etwas wie

distributions = [st.laplace, st.norm, st.expon, st.dweibull, st.invweibull, st.lognorm, st.uniform] 
distributionPairs = [[modelA.name, modelB.name] for modelA in distributions for modelB in distributions] 

und verwenden Sie diese Paare, um einen MLE-Wert dieses Verteilungspaares zu erhalten, der zu den Daten passt?

Answer 1

Es ist keine vollständige Antwort, aber es kann Ihnen helfen, Ihr Problem zu lösen. Sagen wir, du weißt, dass dein Problem durch zwei Dichten erzeugt wird. Eine Lösung wäre K-Mittelwert oder EM-Algorithmus zu verwenden.

Initalisierung. Sie initialisieren Ihren Algorithmus, indem Sie jede Beobachtung auf die eine oder andere Dichte anwenden. Und Sie initialisieren die zwei Dichten (Sie initialisieren die Parameter der Dichte, und einer der Parameter in Ihrem Fall ist "Gaussian", "Laplace", und so weiter ... Iteration. Dann iterately, führen Sie die beiden folgende Schritte:..

Schritt 1. die Parameter optimieren davon aus, dass die Affektiertheit von jedem Punkt richtig ist Sie nun jede Optimierung Solver verwendet Dieser Schritt können Sie mit einer Schätzung der besten zwei Dichten (mit bestimmten Parametern) die zu Ihren Daten passen

Schritt 2. Sie klassifizieren jede Beobachtung zu einem d oder die andere nach der größten Wahrscheinlichkeit.

Sie wiederholen bis zur Konvergenz.

Dies ist sehr gut in dieser Web-Seite erklärt https://people.duke.edu/~ccc14/sta-663/EMAlgorithm.html

Wenn Sie nicht wissen, wie viele Dichten haben Ihre Daten erzeugt, ist das Problem schwieriger. Sie müssen mit dem Problem der Strafe umgehen, das ein bisschen schwieriger ist.

Hier ist ein Coding-Beispiel in einem einfachen Fall: Sie wissen, dass Ihre Daten von 2 verschiedenen Gaussians kommen (Sie wissen nicht, wie viele Variablen aus jeder Dichte generiert werden). In Ihrem Fall können Sie diesen Code in einer Schleife auf jedem möglichen Paar Dichte (rechnerisch länger, aber ich würde empirisch arbeiten vermuten) einstellen

import scipy.stats as st 
import numpy as np 

#hard coded data generation 
data = np.random.normal(-3, 1, size = 1000) 
data[600:] = np.random.normal(loc = 3, scale = 2, size=400) 

#initialization 

mu1 = -1 
sigma1 = 1 

mu2 = 1 
sigma2 = 1 

#criterion to stop iteration 
epsilon = 0.1 
stop = False 

while not stop : 
    #step1 
    classification = np.zeros(len(data)) 
    classification[st.norm.pdf(data, mu1, sigma1) > st.norm.pdf(data, mu2, sigma2)] = 1 

    mu1_old, mu2_old, sigma1_old, sigma2_old = mu1, mu2, sigma1, sigma2 

    #step2 
    pars1 = st.norm.fit(data[classification == 1]) 
    mu1, sigma1 = pars1 
    pars2 = st.norm.fit(data[classification == 0]) 
    mu2, sigma2 = pars2 

    #stopping criterion 
    stop = ((mu1_old - mu1)**2 + (mu2_old - mu2)**2 +(sigma1_old - sigma1)**2 +(sigma2_old - sigma2)**2) < epsilon 

#result  
print("The first density is gaussian :", mu1, sigma1) 
print("The first density is gaussian :", mu2, sigma2) 
print("A rate of ", np.mean(classification), "is classified in the first density") 

Hoffe, es hilft.