Reading Markov Chain Monte Carlo (MCMC) and the numerical integration method to get the posterior.

enter image description here


The implementation code divides the posterior with (post.sum() / len(thetas)). Why the normalising P(X) becomes the (post.sum() / len(thetas))?

thetas = np.linspace(0, 1, 200)
prior = st.beta(a, b)

post = prior.pdf(thetas) * st.binom(n, thetas).pmf(h)
post /= (post.sum() / len(thetas))    # <---------------

Kindly provide intuitive explanation why P(X) is posterior.sum() / (Number of samples).

It seems related with the Monte Carlo Integration where the professor mentioned Dirac Delta property of extracting the spikes in YouTube Machine learning - Logistic regression, but could not understand from where the formula came from.

enter image description here


from __future__ import division
import os
import sys
import glob
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as st

%matplotlib inline
%precision 4

from mpl_toolkits.mplot3d import Axes3D
import scipy.stats as stats
from functools import partial

n = 100
h = 61
p = h/n
rv = st.binom(n, p)
mu = rv.mean()

a, b = 10, 10
prior = st.beta(a, b)
thetas = np.linspace(0, 1, 200)

post = prior.pdf(thetas) * st.binom(n, thetas).pmf(h)
post /= (post.sum() / len(thetas))

plt.figure(figsize=(12, 9))
plt.plot(thetas, prior.pdf(thetas), label='Prior', c='blue')
plt.plot(thetas, n*st.binom(n, thetas).pmf(h), label='Likelihood', c='green')
plt.plot(thetas, post, label='Posterior', c='red')
plt.xlim([0, 1])
plt.xlabel(r'$\theta$', fontsize=14)
plt.ylabel('Density', fontsize=16)

enter image description here



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.