Modelo de ajuste para duas distribuições normais no PyMC

Como sou um engenheiro de software tentando aprender mais estatísticas, você terá que me perdoar antes mesmo de começar, esse é um território novo e sério ...

Estou aprendendo PyMC e trabalhando com alguns exemplos realmente (realmente) simples. Um problema para o qual não consigo trabalhar (e não consigo encontrar exemplos relacionados) é ajustar um modelo aos dados gerados a partir de duas distribuições normais.

Digamos que tenho 1000 valores; 500 gerados de um Normal(mean=100, stddev=20)e outros 500 gerados de um Normal(mean=200, stddev=20).

Se eu quiser ajustar um modelo a eles, ou seja, determine as duas médias e o único desvio padrão, usando PyMC. Eu sei que é algo do tipo ...

mean1 = Uniform('mean1', lower=0.0, upper=200.0)
mean2 = Uniform('mean2', lower=0.0, upper=200.0)
precision = Gamma('precision', alpha=0.1, beta=0.1)

data = read_data_from_file_or_whatever()

@deterministic(plot=False)
def mean(m1=mean1, m2=mean2):
    # but what goes here?

process = Normal('process', mu=mean, tau=precision, value=data, observed=True)

ou seja, o processo de geração é Normal, mas mu é um dos dois valores. Só não sei como representar a "decisão" entre se um valor vem m1ou m2.

Talvez eu esteja apenas adotando completamente a abordagem errada para modelar isso? Alguém pode me indicar um exemplo? Eu posso ler BUGS e JAGS, então tudo está realmente bom.

Você está absolutamente certo de que metade veio de uma distribuição e a outra metade da outra? Caso contrário, podemos modelar a proporção como uma variável aleatória (o que é algo muito bayesiano a se fazer).

A seguir, o que eu faria, algumas dicas são incorporadas.

from pymc import *

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 ) #better to use normals versus Uniforms (unless you are certain the value is  truncated at 0 and 200 
mean2 = Normal( "mean2", 0, 0.001 )

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


#generate some artificial data   
v = np.random.randint( 0, 2, size)
data = v*(10+ np.random.randn(size) ) + (1-v)*(-10 + np.random.randn(size ) )


obs = Normal( "obs", mean, precision, value = data, observed = True)

model = Model( {"p":p, "precision": precision, "mean1": mean1, "mean2":mean2, "obs":obs} )

Alguns pontos, relacionados à discussão acima:

1. A escolha de normal difuso versus uniforme é bastante acadêmica, a menos que (a) você esteja preocupado com a conjugação; nesse caso, você usaria o normal ou (b) há alguma chance razoável de que o verdadeiro valor possa estar fora dos pontos finais do uniforme . Com o PyMC, não há razão para se preocupar com a conjugação, a menos que você deseje especificamente usar um amostrador Gibbs.

2. Na verdade, uma gama não é uma ótima opção para um não informativo antes de um parâmetro de variação / precisão. Pode acabar sendo mais informativo do que você pensa. Uma escolha melhor é colocar um uniforme antes do desvio padrão e depois transformá-lo por um quadrado inverso. Veja Gelman 2006 para detalhes.