0
$\begingroup$

Background

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

enter image description here

Question

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

Code

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
plt.style.use('ggplot')

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

np.random.seed(1234)
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)
plt.legend();

enter image description here

Researches

$\endgroup$

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.