I am new to Bayesian Statistics and pymc3. I follow tutorials to build models, but want to understand more about how pm.sample() do. Let’s use a simple example. The following model is from official tutrial:
basic_model = pm.Model()
with basic_model:
# Priors for unknown model parameters
alpha = pm.Normal('alpha', mu=0, sigma=10)
beta = pm.Normal('beta', mu=0, sigma=10, shape=2)
sigma = pm.HalfNormal('sigma', sigma=1)
# Expected value of outcome
mu = alpha + beta[0]*X1 + beta[1]*X2
# Likelihood (sampling distribution) of observations
Y_obs = pm.Normal('Y_obs', mu=mu, sigma=sigma, observed=Y)
with basic_model:
# draw 500 posterior samples
trace = pm.sample(500)
Is there an explanation how pm.sample(500) work step by step to actually do sampling the posterior distribution defined by basic_model? When and how the priors and likelihood are used in the sampling? What is the difference between pm.Normal() with observed=Y argument and that without such argument? How and when to use them?
As for code logic, you can have a look at the dev guide https://docs.pymc.io/developer_guide.html (but I suggest you to get conformable to Bayesian Concepts first).
Below are some short answers to your questions:
Is there an explanation how pm.sample(500) work step by step to actually do sampling the posterior distribution defined by basic_model?
See https://docs.pymc.io/developer_guide.html#logp-and-dlogp. Basically, pymc3.Model computes everything in the with context and generate a logp_dlogp function, and then a sampler will use that logp_dlogp function to do sampling
When and how the priors and likelihood are used in the sampling?
What is the difference between pm.Normal() with observed=Y argument and that without such argument? How and when to use them?
Without observed=Y it is a free parameter that you want to infer it’s posterior distribution. In principle you use it when you want to associate your data (i.e., observation) to some distribution.