Estimating the true temperature of a room with the Metropolis algorithm

In this notebook I study a simple application of the Metropolis-Hastings algorithm.

Consider a room and a low-cost thermometer. If we measure the temperature several times, the values are not all identical: the sensor introduces a random error.

The true temperature of the room is indicated by $T$ .

The observed measurements are indicated with:

y_1, y_2, y_3, \dots, y_n

The problem is to estimate which values of $T$ are compatible with the collected measurements.

Measurement model

Suppose that each measurement is composed of the true temperature plus a random error:

y_i = T + \varepsilon_i

where $y_i$ is measurement number $i$ , $T$ is the true temperature of the room, and $\varepsilon_i$ is the thermometer error.

We assume that the sensor error is Gaussian:

\varepsilon_i \sim \mathcal{N}(0, \sigma^2)

The parameter $\sigma$ represents the uncertainty of the thermometer.

If $\sigma$ is small, the sensor is precise. If $\sigma$ is large, the measurements fluctuate strongly around the true temperature.

Data generation

To build the simulation I generate a set of artificial measurements.

I choose a true temperature:

T_{\text{true}} = 22.5

and a sensor uncertainty:

\sigma = 0.8

Each measurement is generated as:

y_i = T_{\text{true}} + \varepsilon_i

with:

\varepsilon_i \sim \mathcal{N}(0, \sigma^2)

In the following, I will use the generated measurements as if they were experimental data. The algorithm will not directly use $T_{\text{true}}$ : that value will only be used to check whether the final estimate is reasonable.

Likelihood

Given a possible temperature value $T$ , we can evaluate how compatible it is with the observed measurements.

If the errors are Gaussian, the probability of the observed measurements given $T$ is proportional to:

p(y|T) \propto \exp\left( -\frac{1}{2\sigma^2} \sum_{i=1}^{n}(y_i - T)^2 \right)

Values of $T$ close to the measurements produce a higher probability.

Values of $T$ far from the measurements produce a lower probability.

Distribution to sample

I want to sample the plausible values of the true temperature.

I indicate this distribution as:

p(T|y)

that is, the probability of $T$ after observing the measurements $y$ .

In this simulation I use a very simple initial piece of information: the true temperature must lie in a plausible interval, for example between 15 °C and 30 °C.

Inside this interval, I initially consider all values possible.

Therefore, the distribution to sample is proportional to the likelihood:

p(T|y) \propto p(y|T)

Metropolis algorithm

The algorithm builds a chain of values:

T_0, T_1, T_2, \dots, T_N

Each value in the chain represents a possible true temperature of the room.

At each step, a new value is proposed:

T_{\text{new}} = T_{\text{old}} + \eta

where $\eta$ is a random perturbation:

\eta \sim \mathcal{N}(0, s^2)

The parameter $s$ controls the average size of the proposed steps.

If the new value explains the measurements better, it is accepted. If it explains them worse, it may still be accepted with a certain probability.

Acceptance probability

The probability of accepting the new value is:

\alpha = \min\left( 1, \frac{p(T_{\text{new}}|y)}{p(T_{\text{old}}|y)} \right)

In our case:

p(T|y) \propto p(y|T)

therefore I can use:

\alpha = \min\left( 1, \frac{p(y|T_{\text{new}})}{p(y|T_{\text{old}})} \right)

If the ratio is greater than 1, the new value is always accepted.

If the ratio is less than 1, the new value is accepted with probability equal to the ratio.

Using the logarithm of the probability

To avoid very small numbers, in the code I will use the logarithm of the probability.

The log-likelihood is:

\log p(y|T) = -\frac{1}{2\sigma^2} \sum_{i=1}^{n}(y_i - T)^2

Additive constants independent of $T$ can be ignored, because they cancel out in the acceptance ratio.

The ratio between probabilities becomes:

\frac{p(y|T_{\text{new}})}{p(y|T_{\text{old}})} = \exp\left( \log p(y|T_{\text{new}}) - \log p(y|T_{\text{old}}) \right)

Initial state

At the beginning, the chain may depend strongly on the chosen initial value.

For this reason, the first samples are discarded. After that, the values generated by the chain are used to estimate the distribution of the true temperature.

From these samples I can calculate:

\text{mean}

\text{standard deviation}

\text{interval of plausible values}

What I want to observe

The simulation will produce three main pieces of information.

First I observe the thermometer measurements, which fluctuate around the true temperature.

Then I observe the Metropolis-Hastings chain, that is, the sequence of proposed and accepted temperatures.

Finally I observe the histogram of the samples after burn-in. This histogram represents the distribution of the plausible values of the true temperature.

The sample mean provides an estimate of the room temperature. The width of the distribution indicates the uncertainty of the estimate.

Preparing the numerical method

import numpy as np

def metropolis_hastings(log_pdf, x0, proposal_sigma, n_steps, seed=None):
    rng = np.random.default_rng(seed)
    x = np.asarray(x0, dtype=float)
    samples = np.empty((n_steps,) + x.shape)
    logp = log_pdf(x)
    accepted = 0

    for k in range(n_steps):
        xp = x + proposal_sigma * rng.normal(size=x.shape)
        logp_p = log_pdf(xp)

        if np.log(rng.random()) < logp_p - logp:
            x, logp = xp, logp_p
            accepted += 1

        samples[k] = x

    return samples, accepted / n_steps

Preparing the thermometer data

import matplotlib.pyplot as plt

# Simulated thermometer data

T_true = 22.5
sigma_sensor = 0.8
n_measurements = 30

rng = np.random.default_rng(1)

measurements = T_true + sigma_sensor * rng.normal(size=n_measurements)

print("Measurements:")
print(measurements)

print("Mean of the measurements =", np.mean(measurements))
print("Standard deviation of the measurements =", np.std(measurements))

Measurements:
[22.77646735 23.15729451 22.76434966 21.45747421 23.22428469 22.85709966
 22.07043741 22.96489448 22.79165792 22.735306   22.52273779 22.93737039
 21.91083673 22.36967204 22.11430455 22.97907697 22.53177769 22.2660346
 21.87447323 22.29424621 22.50651374 22.27951768 23.53525105 23.30537945
 20.33107002 20.9887894  22.36018233 22.16224767 22.6709144  22.67385754]
Mean of the measurements = 22.447117312857806
Standard deviation of the measurements = 0.662021471090821

Plot of the measurements

plt.figure(figsize=(9, 5))

plt.plot(measurements, "o", label="Thermometer measurements")
plt.axhline(T_true, linestyle="-", color="red", label="True temperature")
plt.axhline(np.mean(measurements), linestyle=":", color="green", label="Mean of the measurements")

plt.xlabel("Measurement number")
plt.ylabel("Temperature [°C]")
plt.title("Noisy temperature measurements")
plt.grid(True)
plt.legend()
plt.show()

def log_temperature_posterior(T):
    T = np.asarray(T)

    # Plausible interval for the room temperature
    if np.any(T < 15.0) or np.any(T > 30.0):
        return -np.inf

    residuals = measurements - T

    return -0.5 * np.sum((residuals / sigma_sensor)**2)

samples, acceptance_rate = metropolis_hastings(
    log_pdf=log_temperature_posterior,
    x0=20.0,
    proposal_sigma=0.5,
    n_steps=20000,
    seed=2
)

print("Acceptance rate =", acceptance_rate)
print("Sample mean =", np.mean(samples))
print("Sample standard deviation =", np.std(samples))

Acceptance rate = 0.34065
Sample mean = 22.444722105834913
Sample standard deviation = 0.147636597594561

plt.figure(figsize=(9, 5))

plt.plot(samples)

plt.axhline(T_true, linestyle="-", color="red", label="True temperature")
plt.axhline(np.mean(measurements), linestyle=":", color="green", label="Mean of the measurements")

plt.xlabel("Metropolis-Hastings step")
plt.ylabel("Temperature [°C]")
plt.title("Metropolis-Hastings chain")
plt.grid(True)
plt.legend()
plt.show()

Burn-in

burn_in = 2000

posterior_samples = samples[burn_in:]

T_estimate = np.mean(posterior_samples)
T_std = np.std(posterior_samples)

print("Total samples =", len(samples))
print("Samples discarded as burn-in =", burn_in)
print("Samples used =", len(posterior_samples))
print("MH temperature estimate =", T_estimate)
print("MH uncertainty =", T_std)
print("True temperature =", T_true)
print("Mean of the measurements =", np.mean(measurements))

Total samples = 20000
Samples discarded as burn-in = 2000
Samples used = 18000
MH temperature estimate = 22.44513423267025
MH uncertainty = 0.14417158278099254
True temperature = 22.5
Mean of the measurements = 22.447117312857806

plt.figure(figsize=(9, 5))

plt.hist(posterior_samples, bins=40, density=True, edgecolor="black")

plt.axvline(T_true, linestyle="-", color="red", label="True temperature")
plt.axvline(np.mean(measurements), linestyle=":", color="green", label="Mean of the measurements")
plt.axvline(T_estimate, linestyle="--", color="orange", label="MH estimate")

plt.xlabel("Temperature [°C]")
plt.ylabel("Density")
plt.title("Sampled distribution of the true temperature")
plt.grid(True)
plt.legend()
plt.show()

The histogram shows the distribution of temperature values sampled by the Metropolis-Hastings algorithm after burn-in.

The distribution is concentrated around the values most compatible with the thermometer measurements.

The sample mean provides an estimate of the true temperature of the room. The width of the histogram represents the uncertainty of the estimate.

Achtung!

Temperature of a room with Metropolis