La production d'un MLE pour une paire de distributions en python

Question

Ok, donc ma courbe de courant de code de montage comporte une étape qui utilise scipy.stats pour déterminer la distribution à droite sur la base des données,

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

Lorsque des données est un liste de valeurs numériques. Cela fonctionne très bien jusqu'ici pour ajuster des distributions unimodales, confirmées dans un script qui génère aléatoirement des valeurs à partir de distributions aléatoires et utilise curve_fit pour redéfinir les paramètres.

Maintenant, je voudrais rendre le code capable de gérer les distributions bimodale, comme dans l'exemple ci-dessous:

Est-il possible d'obtenir un MLE pour une paire de modèles de scipy.stats afin de déterminer si une paire particulière de distributions est un bon ajustement pour les données ?, quelque chose comme

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] 

et utilisez ces paires pour obtenir une valeur MLE de cette paire de distributions correspondant aux données?

Answer 1

Ce n'est pas une réponse complète mais cela peut vous aider à résoudre votre problème. Disons que vous savez que votre problème est généré par deux densités. Une solution serait d'utiliser k-mean ou l'algorithme EM.

Initalisation. Vous initialisez votre algorithme en affectant chaque observation à l'une ou l'autre densité. Et vous initialisez les deux densités (vous initialisez les paramètres de la densité, et l'un des paramètres dans votre cas est "gaussien", "laplace", et ainsi de suite ... Itération Ensuite, itérativement, vous exécutez les deux procédez comme suit:..

étape 1. optimiser les paramètres en supposant que le vous pouvez maintenant utiliser un solveur d'optimisation afféterie de chaque point est juste Cette étape vous fournir une estimation des deux meilleures densités (avec paramètre donné) qui correspond à vos données

Étape 2. Vous classifiez chaque observation à un d l'ensity ou l'autre selon la plus grande probabilité.

Vous répétez jusqu'à convergence.

Ceci est très bien expliqué dans cette page web https://people.duke.edu/~ccc14/sta-663/EMAlgorithm.html

Si vous ne savez pas combien de densités ont généré vos données, le problème est plus difficile. Vous devez travailler avec un problème de classification pénalisé, ce qui est un peu plus difficile.

Voici un exemple de codage dans un cas facile: vous savez que vos données proviennent de 2 gaussiennes différentes (vous ne savez pas combien de variables sont générées à partir de chaque densité). Dans votre cas, vous pouvez régler ce code à boucle sur toutes les paires possibles de densité (informatiquement plus, mais travailleriez de manière empirique, je présume)

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

Hope it helps.