Full Bayesian Parameter Estimation - sampling

Bayesian parameter estimation

The previous post introduced the basic concept of Bayesian parameter estimation. If you haven't read it, please check that post before this one.

When Bayesian estimation is not simple

There are two cases when Bayesian estimation is tractable. The simplest is when a prior distribution is a conjugate prior. Computing the normalising constant of a single or low-dimensional parameter of \( \theta \) is also tractable with a grid approximation.

In this post, I will stick to the Pikachu encounter rate example like in the previous post too.

Imagine we have analysed 10 routes, each with its own Pikachu encounter rate: \[ \boldsymbol{\theta} = [\theta_1, \dots, \theta_{10}], \quad \theta_k \in [0, 1] \] On each route \( k \), we observed \( n_k \) Pikachu out of \( N_k \) Pokemon, so the likelihood is a product of Binomials, \( P(D|\boldsymbol{\theta}) \propto \prod_{k=1}^{10} \theta_k^{n_k} (1-\theta_k)^{N_k - n_k} \).

We use a single multivariate Gaussian prior over \( \boldsymbol{\theta} \), with mean vector \( \boldsymbol{\mu} \) and covariance matrix \( \boldsymbol{\Sigma} \) that encodes the correlation between similar routes: \[ P(\boldsymbol{\theta}) \propto \exp\left( -\frac{1}{2} (\boldsymbol{\theta} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\boldsymbol{\theta} - \boldsymbol{\mu}) \right) \] This prior is again not conjugate to the Binomial likelihood, so the posterior has no clean solution. The normalisation constant is now a 10-dimensional integral: \[ P(D) = \int_0^1 \cdots \int_0^1 \prod_{k=1}^{10} \theta_k^{n_k} (1-\theta_k)^{N_k - n_k} \exp\left( -\frac{1}{2} (\boldsymbol{\theta} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\boldsymbol{\theta} - \boldsymbol{\mu}) \right) d\theta_1 \cdots d\theta_{10} \] With one parameter, a grid of 1,000 steps meant evaluating the unnormalised posterior 1,000 times. With ten parameters, the grid must cover every combination of values, so the count becomes \( 1000^{10} = 10^{30} \) evaluations. And because \( \boldsymbol{\Sigma} \) correlates the routes, all parameter combinations must be explored together.

This normalising constant with a high dimensional parameter is therefore tough to compute, even with a modern computer.

Sampling

Sampling is a technique to approximate a distribution. As the name suggests, sampling draws samples from a target distribution. The process requires \( \text{likelihood} \times \text{prior} \) without the normalisation constant. The idea is that enough samples form the shape of the posterior distribution. We need neither a conjugate prior nor a grid approximation to estimate a posterior distribution that has no closed form.

Markov Chain Monte Carlo (MCMC)

One type of sampling technique which I'm going to introduce in this post is Markov Chain Monte Carlo (MCMC). "Markov Chain" means that each sample is generated from the previous sample. Samples generated in the past apart from the previous one do not contribute to generation of the current sample. The "Monte Carlo" part refers to the randomness of the method. Essentially, Monte Carlo estimation uses random samples to approximate quantities (integrals, expectations, probabilities, and distributions). It's called Monte Carlo because the casino in Monte Carlo relies on probabilities to steal money gain profit from participants. I think this naming is genius.

Metropolis-Hastings

There are several sampling algorithms that belong to MCMC, but the classic MCMC algorithm I'm going to discuss in this post is the Metropolis-Hastings algorithm. The core concept of the Metropolis-Hastings algorithm is "propose a sample, then accept or reject it".
Could this be the coolest named algorithm?

The goal of the Metropolis-Hastings (MH) algorithm is to generate samples from a distribution \( \pi(\theta) \). \( \pi(\theta) \) is the posterior distribution \( P(\theta|D) \) without the normalising constant.

The MH algorithm randomly chooses an initial parameter \( \theta \) and iterates as follows:
  • Given a \( \theta \), a proposal distribution \( Q(\theta'|\theta) \) proposes the value \( \theta' \) as the next sample.
  • Compute the acceptance probability. \[ A(\theta'|\theta) = \min\left(1, \frac{\pi(\theta')Q(\theta|\theta')}{\pi(\theta)Q(\theta'|\theta)}\right) \] 
    • Accept the proposal with probability of \( A \) and set \( \theta = \theta' \).
    • Reject the proposal with probability of \( 1 - A \) and keep \( \theta = \theta \).
  • Record \( \theta \) as a sample.

Metropolis Hastings: Example

For the demonstration purpose, I'm going to use the posterior distribution of a single dimensional parameter from the previous post: \[ \pi(\theta) = \theta^k \times (1 - \theta)^{n-k} \times \exp\left(-\frac{(\theta-\mu)^2}{2\sigma^2}\right) \] Computing the normalising constant of this posterior distribution is possible with the grid approximation, but let's work out with MH here.

The Binomial likelihood and a truncated Normal distribution are defined as follows in the previous post: 
  • The Binomial likelihood has this data: 1 out of 5 Pokemon is Pikachu \(\binom{5}{1}\theta^1(1-\theta)^4\).
  • The Normal distribution as a prior uses the mean \( \mu = 0.1 \) and the standard deviation \( \sigma = 0.1 \). 
Implementation of HM is something like this. 
n_steps = 50_000
step_size = 0.05        # std-dev of the Gaussian proposal
current = 0.5           # arbitrary starting point
samples = []
for _ in range(n_steps):
    proposal = current + rng.normal(0, step_size)
    ratio = unnormalised_posterior(proposal) / unnormalised_posterior(current)
    if rng.uniform() < ratio:   # uphill: always, downhill: with prob = ratio
        current = proposal
    samples.append(current)     # record even when the proposal is rejected

Observing this code block tells us the following things.
  • The normal distribution "rng.normal(0, step_size)" is the proposal distribution \( Q(\theta'|\theta) \).
  • The acceptance probability is computed only based on unnormalised posteriors \( \frac{\pi(\theta')}{\pi(\theta)} \).
  • "if rng.uniform() < ratio: current = proposal", this is where the current starting point is updated with the proposed value.
Let's look into each element for understanding the HM algorithm.  

    Proposal distribution is the Normal distribution

    The Normal distribution is symmetric and is a convenient choice for the proposal distribution. The Normal distribution here has mean \( \mu=0 \) and standard devication \( \sigma = 0.05 \). Below are three randomly generated example values from this Normal distribution. 
    >>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> rng.normal(0, 0.05)
0.015235853987721568
>>> rng.normal(0, 0.05)
-0.05199920531202478
>>> rng.normal(0, 0.05)
0.03752255970032287
    Visualisation of this distribution looks like the graph below. 
    The proposal distribution \( Q \) in the numerator and the denominator of the acceptance probability \( \frac{Q(\theta|\theta')}{Q(\theta'|\theta)} \) would look like this: \[ Q(\theta|\theta') = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(\theta-\theta')^2}{2\sigma^2} \right) \] \[ Q(\theta'|\theta) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(\theta'-\theta)^2}{2\sigma^2} \right) \] Essentially, the only difference is in the exponent, where \( (\theta - \theta')^2 = (\theta' - \theta)^2 \). This ratio \( \frac{Q(\theta|\theta')}{Q(\theta'|\theta)} \) therefore equals 1.

    In other words, the proposal distribution is symmetric, because the initial value \( \theta \) could take a step of \( \theta' = \theta + 0.01 \) equally likely to a step of \( \theta' = \theta - 0.01 \). We don't need to compute the Q ratio. If the proposal distribution is asymmetric, the full Q ratio must be computed.

    Acceptance probability

    After cancelling the Q ratio, what is left for the acceptance probability is the ratio of posteriors: \[ A(\theta'|\theta) = \min \left(1,  \frac{\pi(\theta')}{\pi(\theta)} \right) \] The unnormalised posterior distribution of the Pikachu encounter rate was defined as follows: \[ \pi(\theta) = \theta^k \times (1 - \theta)^{n-k} \times \exp\left(-\frac{(\theta-\mu)^2}{2\sigma^2}\right) \] Using the concrete values defined above, implementation of the posterior computation is below: 
    n, k = 5, 1
mu, sigma = 0.1, 0.1
a, b = (0 - mu) / sigma, (1 - mu) / sigma

def likelihood(theta):
    return comb(n, k) * theta**k * (1 - theta)**(n - k)

def prior(theta):
    # truncnorm returns 0 outside [0, 1], so out-of-range
    # proposals are rejected automatically
    return truncnorm.pdf(theta, a, b, loc=mu, scale=sigma)

def unnormalised_posterior(theta):
    return likelihood(theta) * prior(theta)

    The ratio \( \frac{\pi(\theta')}{\pi(\theta)} \) compares the posterior density at the proposed parameter \( \theta' \) against the current parameter \( \theta \). If \( \theta' \) has the higher density, the ratio exceeds 1, and since the min caps the acceptance probability at 1, the new parameter \( \theta' \) is accepted. If \( \theta' \) has the lower density, the ratio is less than 1, so \( A \) is equal to the ratio and the proposal is accepted with that probability.

    Value update examples

    Below are the first 5 updates sampled with the MH algorithm.
    proposal: 0.5207738876475073
ratio: 0.37472321449483587
rejected

proposal: 0.5106677452347335
ratio: 0.6080158660811117
accepted: current=0.5106677452347335

proposal: 0.4750292068030035
ratio: 4.9973877117575025
accepted: current=0.4750292068030035

proposal: 0.4325111992605819
ratio: 5.595182105753271
accepted: current=0.4325111992605819

proposal: 0.368528287066306
ratio: 8.93566538454108
accepted: current=0.368528287066306

    The starting value is set to \( \theta = 0.5 \). The highest density of this example posterior distribution appears at \( \theta = 0.131 \).

    The 1st proposal moved \( \theta \) away from the highest density and it was rejected by falling into the 63% of rejection region. The 2nd proposal shifted \( \theta \) away from the highest density again, but this time it was accepted with probability \( 0.61 \). This is when the algorithm explores a proposal, rather than only seeking for values that lead to higher density. From the 3rd to 5th trials, the proposals consistently gave a ratio higher than \(1.0\) and were accepted. 
    Figure.

    Metropolis Hastings: Sampled Distribution

    The full run of the MH algorithm produces samples whose histogram closely matches the the posterior from the grid approximation, as shown below.

    Summary

    This post covered the second topic of Bayesian parameter estimation - when the Bayesian estimation is not tractable. The normalising constant of the posterior distribution can be intractable, when the prior distribution is not conjugate or the parameter dimension is high. Sampling approximates the posterior by drawing samples from the posterior using only the unnormalised form, \( \text{likelihood} \times \text{prior} \). We looked into the Metropolis-Hastings algorithm, which is a member of the Markov chain Monte Carlo (MCMC) algorithm class. 


    Comments

    Post a Comment

    Popular posts from this blog

    Digital Signal Processing Basics

    Image
    Acoustic Signal I had a hard time understanding basics of digital signal processing. I think that the reason was the fact that sounds were not visible. There is a way to "see" the sound, though. This YouTube video for example demonstrates that special cameras make the sound visible.   The sound is the virbation of particles in the air. Something invisible surrounding us rapidly moves back and forth. And our ears can hear this movement. I am writing this post to cover very basics of digital signal processing that I was very slow to understand. I will write concepts of Analogue-to-Digital conversion, sampling and Nyquist frequency. Analogue vs Digital Apart from the sound not usually visible to us, I think that I initially mixed up the concept of the digital signal with the analogue signal. I was still new to the idea that computers represent everything in discrete numbers. Two important things: An anal...
    Read more

    Setup git pre-push hook

    Motivation I often forget to run tests before pushing my commits. This is a bad practice. Someone might use my broken code in a collaborative environment. I was looking for a way to automatically run tests before creating a commit or pushing commits. The solution is Git Hooks:  Customizing Git Hooks . Git Hooks Many hooks can be setup in a Git repository. For example, Git Hooks can run before creating a commit ( pre-commit ), after creating a commit ( post-commit ) and before pushing commits to a remote branch ( pre-push ). All of these are client-side hooks and client-side hooks are not shared across local repositories. The other type of hooks are server-side hooks  that run after commits are pushed to a remote server and policy is imposed on all users of the repository. I have setup  pre-push  in a Git repository to avoid forgetting to run tests before sending updates to the remote server. Setup pre-push The way to setup  pr...
    Read more

    Sound frequency

    Image
    Introduction The signal frequency is the "pitch" of the sound. Some facts about sound frequencies you might encounter in a pub quiz... Typically, male voices range from 85 to 180 Hz and female voices from 165 to 255 Hz .  Humans can easily hear the sound frequency up to 8,000 Hz and lose abilities to hear sounds beyond that frequency through age. The music note C is 261.63 Hz and E 329.63 Hz. The Python notebook is a convenient playground to generate sounds of those frequencies and listen to the sounds. https://github.com/yasumori/blog/blob/main/2025/2025_12_21_signal2.ipynb . An example code snipet to generate a 2,000 Hz sound is also below: import numpy as np from IPython import display def gen_audio(frequency, duration, sample_rate): t = np.linspace(0, duration, duration * sample_rate) return np.sin(2 * np.pi * frequency * t) hz_2000 = gen_audio(2000, 3, 44100) disp...
    Read more