55  πŸ’‘ Estimating Household Energy Savings

A local council introduces a sustainability program encouraging households to reduce electricity use. The program provides households with energy-efficiency advice, LED lighting, and smart-meter feedback.

After one month, the council records the monthly electricity savings (in kWh) for participating households:

\[ Y_i = \text{electricity use before} - \text{electricity use after}. \]

Positive values indicate reduced electricity use.

The council would like to answer the following question:

What is the average monthly electricity saving, and how confident are we that the program reduces electricity use by at least 20 kWh per household?

Bayesian Normal Model

Assume the household savings follow a normal distribution:

\[ Y_i \sim N(\mu,\sigma^2), \qquad i=1,2,\dots,n, \]

where:

  • \(\mu\) is the average household electricity saving;
  • \(\sigma^2\) measures household-to-household variability.

In a Bayesian framework, we must specify prior distributions for the unknown parameters.

A common conjugate prior specification is:

\[ \mu \mid \sigma^2 \sim N(\mu_0,\tau_0^2), \]

\[ \sigma^2 \sim \text{Inverse-Gamma}(\alpha_0,\beta_0). \]

The Normal likelihood together with conjugate priors produces full conditional posterior distributions with standard forms. As a result, Gibbs sampling can be implemented efficiently using direct simulation from familiar distributions.

Full Conditional Posterior Distributions

Using Bayes’ theorem, the resulting full conditional posterior distributions are:

\[ \mu \mid \sigma^2,y \sim N(M_n,V_n), \]

where

\[ V_n= \left( \frac{n}{\sigma^2} + \frac{1}{\tau_0^2} \right)^{-1}, \]

and

\[ M_n= V_n \left( \frac{\sum_{i=1}^n y_i}{\sigma^2} + \frac{\mu_0}{\tau_0^2} \right). \]

Similarly,

\[ \sigma^2 \mid \mu,y \sim \text{Inverse-Gamma}(\alpha_n,\beta_n), \]

where

\[ \alpha_n= \alpha_0+\frac{n}{2}, \]

and

\[ \beta_n= \beta_0 + \frac12 \sum_{i=1}^n (y_i-\mu)^2. \]

These are the conditional posterior distributions used in the Gibbs sampler.

Gibbs Sampling Idea

Instead of sampling directly from the joint posterior distribution of \((\mu,\sigma^2)\), Gibbs sampling repeatedly updates each parameter conditional on the current value of the other.

The Gibbs sampling procedure is:

  1. Choose initial values for \(\mu\) and \(\sigma^2\).
  2. Sample a new value of \(\mu\) given the current value of \(\sigma^2\).
  3. Sample a new value of \(\sigma^2\) given the updated value of \(\mu\).
  4. Repeat many times.

That is,

\[ \mu^{(t+1)} \sim p(\mu \mid \sigma^{2(t)}, y), \]

\[ \sigma^{2(t+1)} \sim p(\sigma^2 \mid \mu^{(t+1)}, y). \]

As the algorithm proceeds, the generated samples approximate draws from the joint posterior distribution of \((\mu,\sigma^2)\).

Gibbs Sampling Implementation

set.seed(1234)

# Simulated household electricity savings (kWh)
y <- c(18, 22, 25, 17, 30, 28, 21, 19, 24, 26,
       15, 23, 27, 20, 29, 31, 16, 22, 24, 18)

n <- length(y)

# Prior hyperparameters
mu0 <- 20        # prior mean
tau2_0 <- 100    # prior variance
alpha0 <- 2      # prior shape
beta0 <- 50      # prior scale

# Gibbs sampler settings
N <- 10000
burnin <- 2000

mu <- numeric(N)
sigma2 <- numeric(N)

# Initial values
mu[1] <- mean(y)
sigma2[1] <- var(y)

for (t in 2:N) {

  # --- Sample mu | sigma2, y ---
  
  var_mu <- 1 / (n / sigma2[t-1] + 1 / tau2_0)

  mean_mu <- var_mu * (
    sum(y) / sigma2[t-1] +
    mu0 / tau2_0
  )

  mu[t] <- rnorm(
    1,
    mean = mean_mu,
    sd = sqrt(var_mu)
  )

  # --- Sample sigma2 | mu, y ---
  
  alpha_post <- alpha0 + n / 2

  beta_post <- beta0 +
    0.5 * sum((y - mu[t])^2)

  sigma2[t] <- 1 / rgamma(
    1,
    shape = alpha_post,
    rate = beta_post
  )
}

# Remove burn-in
mu_post <- mu[(burnin + 1):N]
sigma2_post <- sigma2[(burnin + 1):N]

# Posterior summaries
mean(mu_post)
[1] 22.70689
quantile(mu_post, c(0.025, 0.5, 0.975))
    2.5%      50%    97.5% 
20.46620 22.71309 24.98493 
# Posterior probability of meaningful saving
mean(mu_post > 20)
[1] 0.989875
Python Version 
import numpy as np
np.random.seed(1234)

# Simulated household electricity savings (kWh)
y = np.array([
    18, 22, 25, 17, 30, 28, 21, 19, 24, 26,
    15, 23, 27, 20, 29, 31, 16, 22, 24, 18
])

n = len(y)

# Prior hyperparameters
mu0 = 20        # prior mean
tau2_0 = 100    # prior variance
alpha0 = 2      # prior shape
beta0 = 50      # prior scale

# Gibbs sampler settings
N = 10000
burnin = 2000

mu = np.zeros(N)
sigma2 = np.zeros(N)

# Initial values
mu[0] = np.mean(y)
sigma2[0] = np.var(y, ddof=1)

for t in range(1, N):

    # --- Sample mu | sigma2, y ---

    var_mu = 1 / (n / sigma2[t - 1] + 1 / tau2_0)

    mean_mu = var_mu * (
        np.sum(y) / sigma2[t - 1]
        + mu0 / tau2_0
    )

    mu[t] = np.random.normal(
        loc=mean_mu,
        scale=np.sqrt(var_mu)
    )

    # --- Sample sigma2 | mu, y ---

    alpha_post = alpha0 + n / 2

    beta_post = beta0 + 0.5 * np.sum((y - mu[t])**2)

    # If X ~ Gamma(shape=a, scale=1/rate=b),
    # then 1/X follows an Inverse-Gamma distribution
    sigma2[t] = 1 / np.random.gamma(
        shape=alpha_post,
        scale=1 / beta_post
    )

# Remove burn-in
mu_post = mu[burnin:]
sigma2_post = sigma2[burnin:]

# Posterior summaries
print("Posterior mean of mu:")
print(np.mean(mu_post))
print("\n95% credible interval for mu:")
print(np.percentile(mu_post, [2.5, 97.5]))
print("\nPosterior probability that mu > 20:")
print(np.mean(mu_post > 20))

The posterior mean of \(\mu\) estimates the average household electricity saving, while the credible interval quantifies uncertainty about the true average saving. From the above output, the posterior mean is approximately 22.71 kWh, and the 95% credible interval is approximately (20.47, 24.98) kWh.

The posterior probability, \(P(\mu > 20 \mid y)\), provides a measure of confidence that the program achieves meaningful energy savings. Given the observed data and model assumptions, there is approximately a 98.9875% posterior probability that the program reduces average household electricity use by more than 20 kWh per household per month.

If this probability were low, the council might consider redesigning the program, targeting different households, or introducing stronger incentives.

Trace Plots

Trace plots help assess whether the Markov chain has reached a stable region of the posterior distribution.

par(mfrow = c(2,1), mar = c(3,4,2,1))

plot(mu,
     type = "l",
     ylab = expression(mu),
     main = expression("Trace Plot for " * mu))

plot(sigma2,
     type = "l",
     ylab = expression(sigma^2),
     xlab = "Iteration",
     main = expression("Trace Plot for " * sigma^2))

Python Version 
import matplotlib.pyplot as plt
fig, axs = plt.subplots(2, 1, figsize=(10, 6))
axs[0].plot(mu)
axs[0].set_title("Trace Plot for $\mu$")
axs[0].set_ylabel("$\mu$")
axs[1].plot(sigma2)
axs[1].set_title("Trace Plot for $\sigma^2$")
axs[1].set_xlabel("Iteration")
axs[1].set_ylabel("$\sigma^2$")
plt.tight_layout()
plt.show()

Good trace plots should:

  • fluctuate around a stable region;
  • show no strong trends;
  • explore the posterior distribution well.

Posterior Distribution of \(\mu\)

hist(mu_post,
     breaks = 30,
     probability = TRUE,
     main = expression("Posterior Distribution of " * mu),
     xlab = expression(mu))

abline(v = 20,
       col = "red",
       lwd = 2,
       lty = 2)

Python Version 
plt.hist(mu_post, bins=30, density=True)
plt.axvline(x=20, color='red', linestyle='--', linewidth=2)
plt.title("Posterior Distribution of $\mu$")
plt.xlabel("$\mu$")
plt.ylabel("Density")
plt.show()

The distribution of the posterior samples for \(\mu\) provides insight into the likely values of the average household electricity saving. The histogram shows the shape of the posterior distribution, while the dashed red line indicates the target saving level of 20 kWh. The histogram is approximately normal, reflecting the conjugate normal model and the influence of the data.