Principiante de PyMC: cómo tomar muestras del modelo ajustado

Estoy probando un modelo muy simple: ajustando un Normal donde supongo que sé la precisión, y solo quiero encontrar la media. El siguiente código parece ajustarse a la Normal correctamente. Pero después de ajustar, quiero muestrear del modelo, es decir, generar nuevos datos que sean similares a mi datavariable. Sé que puedo usar trace("mean")para obtener muestras para la variable media. Pero, ¿cómo puedo obtener nuevas muestras del modelo en sí?

He mirado documentos, por ejemplo, http://pymc-devs.github.io/pymc/database.html#accessing-sampled-data . También he visto bastantes ejemplos, por ejemplo, los desastres mineros, y varios de los cuadernos de programación probabilística, y ninguno menciona esto. ¡Yo (más o menos un principiante de MCMC) esperaba que el muestreo del modelo ajustado fuera el punto central! ¿Qué me estoy perdiendo?

from pymc import *
data = np.array([-1, 0, 4, 0, 2, -2, 1, 0, 0, 2, 1, -3, -1, 0, 0, 1, 0, 1])
mean = Uniform("mean", -4, 4)
precision = 2.0**-2
obs = Normal("obs", mean, precision, value=data, observed=True)
model = Model( {"mean": mean, "obs": obs})
mcmc = MCMC(model)
mcmc.sample(10000, 1000, 1)
# I can get samples for the "mean" variable
mean_samples = mcmc.trace("mean")[:]
hist(mean_samples)
# but how can I do the equivalent of mcmc.trace("obs")?

Estás buscando lo que se llama distribución predictiva . Incluir esto es muy simple. Antes de crear el Model, agregue la variable estocástica adicional:

predictive = mc.Normal( "predictive", mean, precision )
model = Model( {"mean": mean, "obs": obs, "pred":predictive})

...

predictive_traces = mcmc.trace("predictive")[:]
hist( predictive_traces )

Esto generará datos artificiales del modelo ajustado. Gracias por llamar mi atención sobre este descuido, lo incluiré en el proyecto BMH.

Llegué aquí varios años después cuando buscaba lo mismo con PyMC3, así que voy a dejar una respuesta relevante para la nueva versión: (de Verificaciones predictivas posteriores ).

ppc = pm.sample_ppc(trace, samples=500, model=model, size=100)

Ahora, ppc contiene 500 conjuntos de datos generados (que contienen 100 muestras cada uno), cada uno con una configuración de parámetros diferente de la posterior.