Solved – Example of “real life” use of Bayesian inference on $\mu$ from a normal distribution

bayesianmeanpriorteaching

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?

Best Answer

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
Related Question