Как смоделировать смесь из 3 нормалей в PyMC?

Question

В CrossValidated возникает вопрос о том, как использовать PyMC для подгонки двух нормальных распределений к данным. Ответ Cosma Colanicchia состоял в том, чтобы использовать распределение Бернулли, чтобы назначить данные одной из двух нормалей:

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-Бернулли. Но как это сделать тогда?

---

редактировать: с предложением CDP, я написал следующий код:

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

E сть mc.Categorical объект, который делает именно это.

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

Возвращает int от 0 до len(p)-1, Для моделирования 3 нормалей, вы делаете p mc.Dirichlet объект (он принимает k длина массива в качестве гиперпараметров; установка значений в массиве одинаковыми - установка предыдущих вероятностей равными). В остальном модель практически идентична.

Это обобщение модели, которую я предложил выше.

---

Обновить:

Итак, вместо того, чтобы использовать разные средства, мы можем свернуть их все в 1:

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

...

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