Bayesian statistics : conjugate prior

from jyquickhelper import add_notebook_menu
add_notebook_menu()

run previous cell, wait for 2 seconds

Reminder

1. Gamma distribution

Solution

1.1. Exercices

We have, in the following a function that generates a gaussian dataset with a gamma prior distribution on the precision. We assume that :

• we know the mean of the gaussian
• we don't know the precision of the gaussian

import numpy as np
import matplotlib.pyplot as plt

def make_normal_dataset(precisions, mu = 5):
    data = np.sqrt(1/precisions) * np.random.randn(len(precisions)) + mu
    return data

# create yout dataset
n = 1000
mu = 5
true_alpha, true_beta = 250, 10

precisions = np.random.gamma(true_alpha, 1/true_beta, size=n)

X = make_normal_dataset(precisions, mu)

# plot the histogram
plt.hist(X, density=True)
plt.xlabel("data values")
plt.ylabel("density")
plt.title(f"histogram of the data")
plt.show()

Our aim is to find out the right precision of the distribution based on the generated dataset X

from scipy import stats

plt.hist(X, density=True)
plt.xlabel("data values")
plt.ylabel("density")
plt.title(f"histogram of the data")
x_range = np.linspace(min(X)-0.5, max(X)+0.5, 100)

# example of precision 2
mu_est, precision_est = 5, 2
pdf_normal = stats.norm.pdf(x_range, mu_est, 1/precision_est)
plt.plot(x_range, pdf_normal, label = f"Normal({mu_est}, {1/precision_est})")

# example of precision 4
mu_est, precision_est = 5, 4
pdf_normal = stats.norm.pdf(x_range, mu, 1/precision_est)
plt.plot(x_range, pdf_normal, label = f"Normal({mu_est}, {1/precision_est})")

plt.legend()
plt.show()

Exercice 1.

As mentioned in the course, in this case the prior and the posterior distributions follow a Gamma distribution. Let's choose $\alpha_{prior} = 1$ and $\beta_{prior} = 1$ as a starting point for our prior. Plot the pdf of this distribution. In the same graph, plot the true distribution.

def gamma_viz(alpha, beta, true_alpha, true_beta):
    # mode
    mode = (alpha-1)/beta
    true_mode = (true_alpha-1)/true_beta

    # plot : estimated
    param_range = np.linspace(0, max(2*mode, 2*true_mode),100)
    pdf_gamma = stats.gamma.pdf(param_range, a=alpha, scale=1/beta)
    plt.plot(param_range, pdf_gamma, label = f"Estimated params : alpha={round(alpha, 1)}, beta={round(beta, 1)}", color="red")
    #plt.axvline(mode, 0, 1, color="red", alpha = 0.3)
    plt.fill_between(param_range, pdf_gamma, color='red', alpha = 0.3)
    
    # plot : true
    param_range = np.linspace(0, max(2*mode, 2*true_mode),100)
    pdf_gamma = stats.gamma.pdf(param_range, a=true_alpha, scale=1/true_beta)
    plt.plot(param_range, pdf_gamma, label = f"True params : alpha={true_alpha}, beta={true_beta}", color="steelblue")
    #plt.axvline(true_mode, 0, 1, color="steelblue", alpha = 0.5)
    plt.fill_between(param_range, pdf_gamma, color='steelblue', alpha = 0.3)
    
    plt.xlabel("parameter values")
    plt.ylabel("density")
    plt.title(f"Distribution of the precision")
    plt.legend()
    return None

• starting point :

alpha_prior = 1
beta_prior = 1
gamma_viz(alpha = alpha_prior, beta = beta_prior, true_alpha = true_alpha, true_beta = true_beta)

Exercice 2. compute numerically $\alpha_{posterior}$ and $\beta_{posterior}$ based on $\alpha_{prior}$, $\beta_{prior}$ and one value of X (the first one for example). plot the estimated distribution along with the true distribution of the precision.

• for one value :

# posterior's parameters
alpha_post = alpha_prior + 1/2
beta_post = beta_prior + np.sum((X[0] - mu)**2)/2

gamma_viz(alpha = alpha_post, beta = beta_post, true_alpha = true_alpha, true_beta = true_beta)

Exercice 3. compute numerically $\alpha_{posterior}$ and $\beta_{posterior}$ based on $\alpha_{prior}$, $\beta_{prior}$ and all values of X. Plot the estimated distribution along with the true distribution of the precision.

• for all values :

# posterior's parameters
alpha_post = alpha_prior + n/2
beta_post = beta_prior + np.sum((X - mu)**2)/2

gamma_viz(alpha = alpha_post, beta = beta_post, true_alpha = true_alpha, true_beta = true_beta)

1.2. Plot's summary

# posterior's parameters
plt.figure(figsize=(16, 10))

for ite, size in enumerate([10, 20, 50, 100, 500, 1000]):
    alpha_post = alpha_prior + size/2
    beta_post = beta_prior + np.sum((X[:size] - mu)**2)/2
    plt.subplot(2, 3, ite+1)
    gamma_viz(alpha = alpha_post, beta = beta_post, true_alpha = true_alpha, true_beta = true_beta)

2. Beta distribution

Solution

3. Gaussian distribution

TODO BEFORE THE END OF THE CLASS : +0.5 point in your final exam if it is correct

4. Advantages

• Easy to do online learning

• We can use the posterior distribution to calculate various quantity/statistics : mean, variance, median, mode, ...

• We can also compute the probability that our parameter lies in an interval

• We can also compute a $T\%$ credible interval for the parameter with $T\in[0, 100]$

5. More conjugate distributions

When the likelihood probability distribution has one of the more common distributions then its conjugate prior can be found in the table of conjugate priors on Wikipedia