BDA3 Chapter 3 Exercise 5

Here’s my solution to exercise 5, chapter 3, of Gelman’s Bayesian Data Analysis (BDA), 3rd edition. There are solutions to some of the exercises on the book’s webpage.

Suppose we weigh an object 5 times with measurements

measurements <- c(10, 10, 12, 11, 9)

all rounded to the nearest kilogram. Assuming the unrounded measurements are normally distributed, we wish to estimate the weight of the object. We will use the uniform non-informative prior \(p(\mu, \log \sigma) \propto 1\).

First, let’s assume the measurments are not rounded. Then the marginal posterior mean is \(\mu \mid y \sim t_{n - 1}(\bar y, s / \sqrt{n}) = t_4(10.4, 0.51)\).

Now, let’s find the posterior assuming rounded measurements. The probability of getting the rounded measurements \(y\) is

\[ p(y \mid \mu, \sigma) = \prod_{i = 1}^n \Phi_{\mu, \sigma} (y_i + 0.5) - \Phi_{\mu, \sigma} (y_i - 0.5) \]

where \(\Phi_{\mu, \sigma}\) is the CDF of the \(\dnorm(\mu, \sigma)\) distribution. This implies that the posterior is

\[ p(\mu, \sigma \mid y) \propto \frac{1}{\sigma^2} \prod_{i = 1}^n \Phi_{\mu, \sigma} (y_i + 0.5) - \Phi_{\mu, \sigma} (y_i - 0.5) . \]

Calculating this marginal posterior mean is pretty difficult, so we’ll use Stan to draw samples. My first attempt at writing the model was a direct translation of the maths above. However, it doesn’t allow us to infer the unrounded values, as required in part d. The model can be expressed differently by considering the unrounded values as uniformly distributed around the rounded values, i.e. \(z_i \sim \duniform (y_i - 0.5, y_i + 0.5)\).

model <- rstan::stan_model('src/ex_03_05_d.stan')

model

S4 class stanmodel 'ex_03_05_d' coded as follows:
data {
  int<lower = 1> n;
  vector[n] y; // rounded measurements
}

parameters {
  real mu; // 'true' weight of the object
  real<lower = 0> sigma; // measurement error
  vector<lower = -0.5, upper = 0.5>[n] err; // rounding error
}

transformed parameters {
  // unrounded values are the rounded values plus some rounding error
  vector[n] z = y + err; // unrounded measurements
}

model {
  target += -2 * log(sigma); // prior
  z ~ normal(mu, sigma);
  // other parameters are uniform
} 

Note that Stan assumes parameters are uniform on their range unless specified otherwise.

Let’s also load a model that assumes the measurements are unrounded.

model_unrounded <- rstan::stan_model('src/ex_03_05_unrounded.stan')

model_unrounded

S4 class stanmodel 'ex_03_05_unrounded' coded as follows:
data {
  int<lower = 1> n;
  vector[n] y; 
}

parameters {
  real mu; 
  real<lower = 0> sigma; 
}

model {
  target += -2 * log(sigma); 
  y ~ normal(mu, sigma);
} 

Now we can fit the models to the data.

data  = list(
  n = length(measurements),
  y = measurements
)
 
fit <- model %>% 
  rstan::sampling(
    data = data,
    warmup = 1000,
    iter = 5000
  ) 

fit_unrounded <- model_unrounded %>% 
  rstan::sampling(
    data = data,
    warmup = 1000,
    iter = 5000
  ) 

We’ll also need some draws from the posteriors to make our comparisons.

draws <- fit %>% 
  tidybayes::spread_draws(mu, sigma, z[index]) %>% 
  # spread out z's so that
  # there is one row per draw
  ungroup() %>%  
  mutate(
    index = paste0('z', as.character(index)),
    model = 'rounded'
  ) %>% 
  spread(index, z)

draws_unrounded <- fit_unrounded %>% 
  tidybayes::spread_draws(mu, sigma) %>% 
  mutate(model = 'unrounded') 

draws_all <- draws %>% 
  bind_rows(draws_unrounded)

The contour plots look very similar but with \(\sigma\) shifted upward when we treat the observations as unrounded measurements. This is contrary to my intuition about what should happen: by introducing uncertainty into our measurments, I would have thought we’d see more uncertainty in our parameter estimates.

The density for \(\mu \mid y\) look much the same in both models. This is expected because the rounded measurement is the mean of all possible unrounded measurements.

The marginal posterior for \(\sigma\) again shows a decrease when taking rounding error into account. I’m not sure why that would happen.

Finally, let’s calculate the posterior for \(\theta := (z_1 - z_2)^2\) (assuming we observe rounded measurements).

sims <- draws %>% 
  mutate(theta = (z1 - z2)^2) 

There is a lot of mass near 0 because the observed rounded measurments are the same for \(z_1\) and \(z_2\). The probability density is also entirely less than 1 because the rounding is off by at most 0.5 in any direction.