CamDavidsonPilon / Probabilistic-Programming-and-Bayesian-Methods-for-Hackers

aka "Bayesian Methods for Hackers": An introduction to Bayesian methods + probabilistic programming with a computation/understanding-first, mathematics-second point of view. All in pure Python ;)
http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/
MIT License
26.66k stars 7.86k forks source link

Ch2 Understanding difference between sampling random() and mcmc.sample() #247

Open rriv opened 9 years ago

rriv commented 9 years ago

Hello,

first thanks a lot to Cameron to bring us such a nice book on MCMC.

Can anyone please tell me if I'm understanding correctly about Parent and Child relationships and distribution sampling.

Let me do the following (from Ch2 example) :

parameter = pm.Exponential("poisson_param", 1)
data_generator = pm.Poisson("data_generator", parameter)

if I do now :

samples1 = [data_generator.random() for i in range(2000)]

I'm taking 2000 random values from data_generator Poisson distribution, but with a constant parameter (initially random chosen from Exponential distribution, but the same for all 2000 sample). I'm telling this following this sentence from PyMC Documentation § 3.4.1 :

Whenever a variable is used as a parent for a child variable, PyMC replaces it with its value attribute when the child’s value or log-probability is computed.

After all, samples1 follows a Poisson distribution. With some parameter, but a Poisson distribution.

If I do now :

model = pm.Model([data_generator,parameter]);
mcmc = pm.MCMC(model)
mcmc.sample(2000);
samples2 = mcmc.trace('data_generator')[:];

Here parameter is drawn each time data_generator is sampled. Now samples2 has no more anything to do with a Poisson distribution. We really mixed the parent random parameter with the child distribution sampling.

So I'd be glad to hear if I'm right or wrong... Thanks in advance,

Robert

CamDavidsonPilon commented 9 years ago

You are correct: the former is just Poisson with a fixed parameter, the latter is a compound distribution of an Exponential and a Poisson.

Perhaps the best way to see this is to plot the samples (using your code, and calling hist on them)

bmh

We can see the differences in the distributions quite clearly (ignoring the width, but focusing on the height of each bar).