I'm working through a book called Bayesian Analysis in Python. The book focuses heavily on the package PyMC3 but is a little vague on the theory behind it.
Say I'm looking at a model like this
My understanding is as such:
To fit the model means to find parameters for "Y~ Normal(mu,sigma)" aka to find good estimates for parameters mu and sigma.
To find each the algorithm samples (randomly) from the prior distributions (for mu that's Normal(M,S), sigma it's HalfCauchy(G))
The metropolis algorithm then decides whether or not to accept it based on (and here is where my understanding gets fuzzy) some likelihood ratio, which i am kinda taking as being the likelihood of the sampled value showing up in distribution y (aka Normal(mu,sigma).pdf(sampled_value)
). This however doesn't make sense, as the probability of mu showing up in that distribution is determined by mu itself... so clearly i'm missing something.
Can someone clarify what is going on in this step?
Questions:
- Am i correct in saying that the algorithm sample randomly from the prior distribution and not from the observed data itself?
- How does the algorithm decide whether to accept or reject the sample?
Bonus points for specific example and/or pseudo code.
Here is the example given by the book for reference:
data = np.array([51.06, 55.12, 53.73, 50.24, 52.05, 56.40, 48.45, 52.34, 55.65, 51.49, 51.86, 63.43, 53.00, 56.09, 51.93, 52.31, 52.33, 57.48, 57.44, 55.14, 53.93, 54.62, 56.09, 68.58, 51.36, 55.47, 50.73, 51.94, 54.95, 50.39, 52.91, 51.5, 52.68, 47.72, 49.73, 51.82, 54.99, 52.84, 53.19, 54.52, 51.46, 53.73, 51.61, 49.81, 52.42, 54.3, 53.84, 53.16])
with pm.Model() as model_g:
mu = pm.Normal('mu', mu=53, sd=2)
sigma = pm.HalfCauchy('sigma', beta=10)
y = pm.Normal('y', mu=mu, sd=sigma, observed=data)
trace_g = pm.sample(1100)