Example of "real life" use of Bayesian inference on $\mu$ from a normal distribution?

Question

A classic example for students, when teaching Bayesian statistics, is to make inference on the mean parameter $\mu$ of a normal distribution, when it has a prior normal distribution.

I would like to find a real-life example of why such a construction would be not only elegant, but also useful.

Can anyone direct me to such an example?

Answer 1

From my perspective the Normal distribution is useful as a prior on the the mean of a Normal distribution because:

1. The Normal is the conjugate prior to the mean of a Normal distribution. Using conjugate priors can really speed up computation in, for example, JAGS

2. It is pretty easy to make the Normal distribution both informative, the mean and SD are parameters that are easy to set as it is relativley clear how they affect the shape of the distribution, and non-informative, just make the SD extremely large and the prior will approach the uniform distribution.

Here is a silly but "real life" example of how the Normal could be used as a slightly informative distribution. :)

Farmer Jöns and Milk Production.

Farmer Jöns has a huge number of cows. He would like to know how much milk a cow is expected to produce on average and thus measures the number of liters of milk that six cows produce during one month:

milk <- c(651, 374, 601, 401, 767, 709)

He sends this data to his statistician friend who tells him that, "Well it's not that much data to work with, but the estimate might get a bit better if you tell me how much milk in your experience a cow produces per month on average". "I haven't thought that much about it", replies Jöns, "but perhaps around 500-600 liters/month is usual."

The statistician decides to run a Bayesian analysis where milk is assumed to be normally distributed:

$$ \mathrm{milk} \sim \mathrm{Norm}(\mu, \sigma) $$

The statistician encodes the information he got from Jöns in the prior on $\mu$ which is pretty easy to do as he choose to use a Normal distribution for this:

$$ \mu \sim \mathrm{Norm}(550, 100) $$

This Normal prior is centered around 550 and has a standard deviation of 100 basically saying that prior to looking at the data the mean is most likely to be in the range 450-650.

As he forgot to ask about the spread of milk production of different cows he puts a flat prior on $\sigma$:

$$ \sigma \sim \mathrm{Unif}(0, 1000) $$

An implementation of this analysis in R using the rjags package would look like this:

library(rjags)

model_string <- "model {
  for(i in 1:length(milk)) {
    milk[i] ~ dnorm(mu, 1/(sigma * sigma) )
  }
  mu ~ dnorm(550, 1 / (100*100))
  sigma ~ dunif(0, 1000)\n
}"
# The reason for the 1/ (sigma*sigma) 'thing' is because dnorm in JAGS is
# parameterized by precision (1 / SD^2) instead of SD as in R.

model <- jags.model(textConnection(model_string), data = list(milk = milk))
fit <- coda.samples(model, variable.names = c("mu", "sigma"), n.iter = 30000)
summary(fit)

 

## 
## Iterations = 1001:31000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 30000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean   SD Naive SE Time-series SE
## mu     570 63.3    0.365          0.368
## sigma  208 89.3    0.516          1.173
## 
## 2. Quantiles for each variable:
## 
##       2.5% 25% 50% 75% 97.5%
## mu     440 530 571 611   693
## sigma  106 150 187 241   435

Looking at the quantiles for the posterior of $\mu$ the statistician tells Jöns that a good guess for the number of liters of milk a cow is expected to produce on average is 571 liters/month and the average is probably not lower than 440 liters/month or higher than 693 liters/month.

What we gained from using a normal distribution as the prior on $\mu$ was a slightly tighter estimate (and hopefully better too). Compare with the 95 % CI using t.test:

t.test(milk)

 

## 
##  One Sample t-test
## 
## data:  milk
## t = 8.819, df = 5, p-value = 0.0003113
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  413.7 754.0
## sample estimates:
## mean of x 
##     583.8

Answer 2

Before delving directly into an example, I'd like to review some of the math for Normal-Normal Bayesian data models. Consider a random sample of n continuous values denoted by $y_1, ..., y_n$. Here the vector $y = (y_1, ..., y_n)^T$ represents the data gathered. The probability model for Normal data with known variance and independent and identically distributed (i.i.d.) samples is

$$ y_1, ..., y_n | \theta \sim N(\theta, \sigma^2) $$

Or as more typically written by Bayesian,

$$ y_1, ..., y_n | \theta \sim N(\theta, \tau) $$

where $\tau = 1 / \sigma^2$; $\tau$ is known as the precision

With this notation, the density for $y_i$ is then

$$ f(y_i | \theta, \tau) = \sqrt(\frac{\tau}{2 \pi}) \times exp\left( -\tau (y_i - \theta)^2 / 2 \right) $$

Classical statistics (i.e. maximum likelihood) gives us an estimate of $\hat{\theta} = \bar{y}$

In a Bayesian perspective, we append maximum likelihood with prior information. A choice of priors for this Normal data model is another Normal distribution for $\theta$. The Normal distribution is conjugate to the Normal distribution.

$$ \theta \sim N(a,1/b) $$

The posterior distribution we obtain from this Normal-Normal (after a lot of algebra) data model is another Normal distribution (so many Normals!)

$$ \theta | y \sim N(\frac{b}{b + n\tau} a + \frac{n \tau}{b + n \tau} \bar{y}, \frac{1}{b + n\tau}) $$

The posterior precision is $b + n\tau$ and mean is a weighted mean between $a$ and $\bar{y}$, $\frac{b}{b + n\tau} a + \frac{n \tau}{b + n \tau} \bar{y}$.

To your comment

I would like to find a real-life example of why such a construction would be not only elegant, but also useful.

The usefulness comes from the fact that you obtain a distribution of $\theta | y$ rather than just an estimate since $\theta$ is viewed as a random variable rather than a fixed (unknown) value. In addition, your estimate of $\theta$ in this model is a weighted average between the empirical mean and prior information.

That said, you can now use any Normal-data textbook example to illustrate this. I'll use the data set airquality within R. Consider the problem of estimating average wind speeds (MPH).

> ## New York Air Quality Measurements
> 
> help("airquality")
> 
> ## Estimating average wind speeds
> 
> wind = airquality$Wind
> hist(wind, col = "gray", border = "white", xlab = "Wind Speed (MPH)")
>

> n = length(wind)
> ybar = mean(wind)
> ybar
[1] 9.957516 ## "frequentist" estimate
> tau = 1/sd(wind)
> 
> 
> ## but based on some research, you felt avgerage wind speeds were closer to 12 mph
> ## but probably no greater than 15,
> ## then a potential prior would be N(12, 2)
> 
> a = 12
> b = 2
> 
> ## Your posterior would be N((1/))
> 
> postmean = 1/(1 + n*tau) * a + n*tau/(1 + n*tau) * ybar
> postsd = 1/(1 + n*tau)
> 
> set.seed(123)
> posterior_sample = rnorm(n = 10000, mean = postmean, sd = postsd)
> hist(posterior_sample, col = "gray", border = "white", xlab = "Wind Speed (MPH)")
> abline(v = median(posterior_sample))
> abline(v = ybar, lty = 3)
> 

> median(posterior_sample)
[1] 10.00324
> quantile(x = posterior_sample, probs = c(0.025, 0.975)) ## confidence intervals
2.5%     97.5% 
9.958984 10.047404 

In this analysis, the researcher (you) can say that given data + prior information, your estimate of average wind, using the 50th percentile, speeds should be 10.00324, greater than simply using the average from the data. You also obtain a full distribution, from which you can extract a 95% credible interval using the 2.5 and 97.5 quantiles.

Below I include two references, I highly recommend reading Casella's short paper. It's specifically aimed at empirical Bayes methods, but explains the general Bayesian methodology for Normal models.

References:

1. Casella, G. (1985). An Introduction to Empirical Bayes Data Analysis. The American Statistician, 39(2), 83-87.

2. Gelman, A. (2004). Bayesian data analysis (2nd ed., Texts in statistical science). Boca Raton, Fla.: Chapman & Hall/CRC.