पीईएमसी में 3 सामान्य के मिश्रण का मॉडल कैसे करें?

Question

डेटा पर दो सामान्य वितरण फिट करने के लिए पीएमएमसी का उपयोग करने के तरीके पर a question on CrossValidated है। Cam.Davidson.Pilon का जवाब दो नॉर्मल्स से एक के लिए डेटा आवंटित करने के लिए एक Bernoulli वितरण का इस्तेमाल किया गया:

size = 10 
p = Uniform("p", 0 , 1) #this is the fraction that come from mean1 vs mean2 
ber = Bernoulli("ber", p = p, size = size) # produces 1 with proportion p. 
precision = Gamma('precision', alpha=0.1, beta=0.1) 

mean1 = Normal("mean1", 0, 0.001) 
mean2 = Normal("mean2", 0, 0.001) 

@deterministic 
def mean(ber = ber, mean1 = mean1, mean2 = mean2): 
    return ber*mean1 + (1-ber)*mean2 

अब मेरी सवाल यह है: कैसे साथ तीन नॉर्मल्स यह करने के लिए?

असल में, मुद्दा यह है कि आप अब बर्नौली वितरण और 1-बर्नौली का उपयोग नहीं कर सकते हैं। लेकिन फिर यह कैसे करें?

---

संपादित करें:

import numpy as np 
import pymc as mc 

n = 3 
ndata = 500 

dd = mc.Dirichlet('dd', theta=(1,)*n) 
category = mc.Categorical('category', p=dd, size=ndata) 

precs = mc.Gamma('precs', alpha=0.1, beta=0.1, size=n) 
means = mc.Normal('means', 0, 0.001, size=n) 

@mc.deterministic 
def mean(category=category, means=means): 
    return means[category] 

@mc.deterministic 
def prec(category=category, precs=precs): 
    return precs[category] 

v = np.random.randint(0, n, ndata) 
data = (v==0)*(50+ np.random.randn(ndata)) \ 
     + (v==1)*(-50 + np.random.randn(ndata)) \ 
     + (v==2)*np.random.randn(ndata) 
obs = mc.Normal('obs', mean, prec, value=data, observed = True) 

model = mc.Model({'dd': dd, 
       'category': category, 
       'precs': precs, 
       'means': means, 
       'obs': obs}) 

निम्नलिखित नमूना प्रक्रिया के साथ निशान के साथ-साथ अच्छे लग रहे हैं: सीडीपी के सुझाव के साथ, मैं निम्नलिखित कोड लिखा था। हल किया!

mcmc = mc.MCMC(model) 
mcmc.sample(50000,0) 
mcmc.trace('means').gettrace()[-1,:] 

Answer 1

mc.Categorical ऑब्जेक्ट है जो यह करता है।

p = [0.2, 0.3, .5] 
t = mc.Categorical('test', p) 
t.random() 
#array(2, dtype=int32) 

यह 0 और len(p)-1 के बीच एक int देता है। 3 सामान्य मॉडल बनाने के लिए, आप pmc.Dirichlet ऑब्जेक्ट बनाते हैं (यह k लंबाई सरणी को हाइपरपरैमीटर के रूप में स्वीकार करता है; सरणी में मानों को सेट करना समान संभावनाओं को समान होने के लिए सेट करना है)। शेष मॉडल लगभग समान है।

यह ऊपर दिए गए मॉडल का एक सामान्यीकरण है।

---

अद्यतन:

ठीक है, तो बजाय विभिन्न साधनों होने के, हम उन सभी 1 टुकड़ों में टूट कर सकते हैं:

means = Normal("means", 0, 0.001, size=3) 

... 

@mc.deterministic 
def mean(categorical=categorical, means = means): 
    return means[categorical]