Metropolitan Sample

Question

I am working on a book called Bayesian Analysis in Python . The book focuses mainly on the PyMC3 package, but is a bit vague in the theory behind it, so I'm pretty confused.

Let's say I have this data:

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

And I look at a model like this:

Using a metropolitan sample , how can I put a model that evaluates mu and sigma.

Here is my hunch about the pseudo code from what I read:

M, S = 50, 1
G = 1

# These are priors right?
mu = stats.norm(loc=M, scale=S)
sigma = stats.halfnorm(scale=G)

target = stats.norm

steps = 1000

mu_samples = [50]
sigma_samples = [1]

for i in range(steps):
    # proposed sample...
    mu_i, sigma_i = mu.rvs(), sigma.rvs()

    # Something happens here
    # How do I calculate the likelidhood??
    "..."
    # some evaluation of a likelihood ratio??

    a = "some"/"ratio"

    acceptance_bar = np.random.random()

    if a > acceptance_bar:
        mu_samples.append(mu_i)
        sigma_samples.append(sigma_i)

What am I missing?

Answer 1

Hope the following example helps you. In this example, I assume that we know the value sigma, so we only have a preliminary option for mu.

sigma = data.std() # we are assuming we know sigma

steps = 1000
mu_old = data.mean() # initial value, just a good guest
mu_samples = []

# we evaluate the prior for the initial point
prior_old = stats.norm(M, S).pdf(mu_old)
# we evaluate the likelihood for the initial point
likelihood_old = np.prod(stats.norm(mu_old, sigma).pdf(data))
# Bayes' theorem (omitting the denominator) for the initial point
post_old = prior_old * likelihood_old

for i in range(steps):
    # proposal distribution, propose a new value from the old one
    mu_new = stats.norm.rvs(mu_old, 0.1)

    # we evaluate the prior
    prior_new = stats.norm(M, S).pdf(mu_new)

    # we evaluate the likelihood
    likelihood_new = np.prod(stats.norm(mu_new, sigma).pdf(data))

    # Bayes' theorem (omitting the denominator)
    post_new = prior_new * likelihood_new

    # the ratio of posteriors (we do not need to know the normalizing constant)
    a =  post_new / post_old

    if np.random.random() < a:
        mu_old = mu_new
        post_old = post_new

    mu_samples.append(mu_old)

Notes:

• , , mu_old 0,1 ( ). MH , PyMC3 ( MH) .
• pdf , logpdf. .
• , , -
• , () .

. , PyMC3 pm.sample() PyMC3 .