Ajuste del modelo para dos distribuciones normales en PyMC

Question

Como soy un ingeniero de software que intenta aprender más sobre estadísticas, tendrás que perdonarme antes de que empiece, esto es territorio de novatos...

He estado aprendiendo PyMC y trabajando con algunos ejemplos muy (muy) sencillos. Un problema que no consigo que funcione (y para el que no encuentro ningún ejemplo relacionado) es el ajuste de un modelo a datos generados a partir de dos distribuciones normales.

Digamos que tengo 1000 valores; 500 generados a partir de un Normal(mean=100, stddev=20) y otros 500 generados a partir de un Normal(mean=200, stddev=20) .

Si quiero ajustar un modelo a ellos, es decir, determinar las dos medias y la desviación estándar única, utilizando PyMC. Sé que es algo en la línea de ...

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)

es decir, el proceso de generación es Normal, pero mu es uno de los dos valores. Sólo que no sé cómo representar la "decisión" entre si un valor proviene de m1 o m2 .

¿Quizás estoy adoptando un enfoque completamente equivocado para modelar esto? ¿Puede alguien indicarme un ejemplo? Puedo leer BUGS y JAGS así que cualquier cosa está bien realmente.

Answer 1

Un par de puntos, relacionados con la discusión anterior:

1. La elección de la normal difusa frente a la uniforme es bastante académica a no ser que (a) te preocupe la conjugación, en cuyo caso utilizarías la normal o (b) haya alguna posibilidad razonable de que el valor verdadero pueda estar fuera de los puntos extremos de la uniforme. Con PyMC, no hay razón para preocuparse por la conjugación, a menos que quiera utilizar específicamente un muestreador de Gibbs.

2. En realidad, una gamma no es una gran elección para un previo no informativo de un parámetro de varianza/precisión. Puede acabar siendo más informativa de lo que se piensa. Una mejor elección es poner una prioridad uniforme en la desviación estándar, y luego transformarla por un cuadrado inverso. Véase Gelman 2006 para más detalles.

Answer 2

¿Está absolutamente seguro de que la mitad proviene de una distribución y la otra mitad de la otra? Si no es así, podemos modelar la proporción como una variable aleatoria (algo muy bayesiano).

Lo siguiente es lo que yo haría, algunos consejos están incrustados.

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