अजगर में वितरण की एक जोड़ी के लिए एक MLE उत्पादन

Question

ठीक है, इसलिए मेरे वर्तमान वक्र ढाले कोड एक कदम scipy.stats का उपयोग करता है सही आंकड़ों के आधार पर वितरण निर्धारित करने के लिए,

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

डेटा कहाँ एक है संख्यात्मक मूल्यों की सूची। यह असामान्य वितरण को फिट करने के लिए अब तक बहुत अच्छा काम कर रहा है, एक ऐसी स्क्रिप्ट में पुष्टि हुई है जो यादृच्छिक रूप से यादृच्छिक वितरण से मूल्य उत्पन्न करती है और पैरामीटर को फिर से निर्धारित करने के लिए curve_fit का उपयोग करती है।

अब मैं कोड bimodal वितरण को संभालने में सक्षम बनाने के लिए चाहते हैं, नीचे दिए गए उदाहरण की तरह:

इसे से मॉडल की एक जोड़ी के लिए एक MLE प्राप्त करने के लिए संभव है scipy.stats यह निर्धारित करने के लिए कि क्या वितरण की एक विशेष जोड़ी डेटा के लिए उपयुक्त है ?,

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] 

और उन जोड़ों का उपयोग डेटा को फ़िट करने वाले वितरण की उस जोड़ी का एमएलई मान प्राप्त करने के लिए करें?

Answer 1

यह एक पूरा उत्तर नहीं है लेकिन यह आपकी समस्या का समाधान करने में आपकी मदद कर सकता है। मान लें कि आप जानते हैं कि आपकी समस्या दो घनत्व से उत्पन्न होती है। एक समाधान के-माध्य या ईएम एल्गोरिदम का उपयोग करना होगा।

इंटाइलाइजेशन। आप प्रत्येक अवलोकन को एक या दूसरे घनत्व को प्रभावित करके अपना एल्गोरिदम प्रारंभ करते हैं। और आप दो घनत्व शुरू करते हैं (आप घनत्व के मानकों को प्रारंभ करते हैं, और आपके मामले में पैरामीटर में से एक "गाऊशियन", "लेपलेस" है, और इसी तरह ... इटरेशन फिर, तुरंत, आप दोनों को चलाते हैं निम्न चरणों का पालन:।।

चरण 1. यह सोचते हैं कि हर बिंदु के दिखावा सही है मापदंडों का अनुकूलन अब आप किसी भी अनुकूलन solver उपयोग कर सकते हैं यह कदम (सबसे अच्छा दो घनत्व की एक अनुमान के साथ आप प्रदान दिया पैरामीटर के साथ) जो आपके डेटा को फिट करता है।

चरण 2. आप प्रत्येक अवलोकन को एक डी में वर्गीकृत करते हैं सबसे बड़ी संभावना के अनुसार उपस्थिति या दूसरा।

आप अभिसरण तक दोहराते हैं।

यह बहुत अच्छी तरह से इस वेब-पेज https://people.duke.edu/~ccc14/sta-663/EMAlgorithm.html

आप नहीं जानते कि कितने घनत्व अपने डेटा जेनरेट किया है से समझाया गया है, समस्या अधिक कठिन है। आपको दंडित वर्गीकरण समस्या के साथ काम करना है, जो थोड़ा कठिन है।

यहां एक आसान मामले में एक कोडिंग उदाहरण है: आप जानते हैं कि आपका डेटा 2 अलग-अलग गॉसियनों से आता है (आप नहीं जानते कि प्रत्येक घनत्व से कितने चर उत्पन्न होते हैं)। आपके मामले में, आप पाश के लिए इस कोड घनत्व के हर संभव जोड़ी पर (computationally अब, लेकिन अनुभव मुझे लगता है काम करेगा) समायोजित कर सकते हैं

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

आशा है कि यह मदद करता है।