Demonstrating numerical accuracy of MCMC on local machines using standard reference data

Modeling
, ,
1

I am trying to show the numerical accuracy of MCMC estimation on my machine. To accomplish this I am using the NIST Statistical Reference Datasets. I am using data from MCMC example 6, which is a very simple model of the form:

y_{i} = mu + epsilon_{i}
epsilon_{i} ~ N(0, sigma^2)

After fitting the model in stan I want to compare to the certified values.

Example code is here:


library(rstan)

response <- read.delim(
             file       = paste0('PATH/mcmc06.txt'),
             skip       = 60,
             sep        = '',
             colClasses = 'numeric',
             header     = FALSE,
             col.names  = 'y'
             )

mcmc6_mod <-
  "
  data {
        int<lower = 0> N;
        real y[N];
        }
        
  parameters {
              real mu;
              real epsilon;
              real<lower = 0> sigma;
              }
              
  model {
         epsilon ~ normal(0, sigma);
         y ~ normal(mu, epsilon);
         }
  "

mcmc6dat <- 
  list(N = nrow(response),y = response$y)

mcmc6_est <-
  rstan::stan(
       model_code = mcmc6_mod,
       model_name = 'NIST MCMC Example6 ',
       data       = mcmc6dat,
       iter       = 10000,
       chains     = 5,
       verbose    = FALSE
       )

print(mcmc6_est )

I’m getting several errors estimating this model:

  • Divergent transitions after warmup
  • Transitions after warmup that exceed max treedepth
  • Large R-hat
  • Bulk ESS too low
  • Tail ESS too low

So, I increase iterations and chains and get similar errors and summaries of parameter estimates that aren’t close.

Then I try simplifying the model, which isn’t ultimately going to work for the purpose of assessing the technical accuracy of my machine with the pre-specified NIST model. I get closer to matching the certified values., but still getting a lot of errors:

mcmc6_mod <-
  "
  data {
        int<lower = 0> N;
        real y[N];
        }
        
  parameters {
              real mu;
              real<lower = 0> epsilon;
              }
              
  model {
         y ~ normal(mu, epsilon);
         }
  "

I get close to matching the moments and quintiles for the estimate mean, but the estimated SD is still way off and I’m getting similar errors.

I’ve attached the datafile here. mcmc06.txt (1.9 KB)

Thanks for any help getting me unstuck.

1 Like
2

I believe the issue is that your Stan model and your mathematical models are not consistent.

Your mathematical model is:

y_{i} = \mu + \epsilon_{i}
\epsilon_{i} \sim N(0, \sigma^2)

But your Stan model is:

y_{i} \sim N(\mu , \epsilon_{i}^2)
\epsilon_{i} \sim N(0, \sigma^2)

The Stan model that you want is:

data {
  int<lower = 0> N;
  real y[N];
 }
        
parameters {
  real mu;
  real<lower = 0> sigma;
}
              
model {
  y ~ normal(mu, sigma);
}
2 Likes
3

The main clarification here is that:

y_{i} = \mu + \epsilon_{i}
\epsilon_{i} \sim N(0, \sigma^2)

Is the same as:

y_{i} \sim N(\mu, \sigma^2)
2 Likes
4

Thank you for clearing up my misunderstanding of the implied model. It’s simpler than I thought.

I did try and run the model in your code before posting, I’ve updated my post with what I tried. I was still getting warnings, but I think given the data I am trying to use here this might just be an issue of increasing chains and iterations, as I tried that this AM and am getting closer with one error that I don’t believe is critical to estimating the model:

Warning message:
In system(paste(CXX, ARGS), ignore.stdout = TRUE, ignore.stderr = TRUE) :
  '-E' not found
5

That warning is safe to ignore, and should be fixed in the next release.

I would guess that the numerical differences are now likely coming from your priors. With only 11 observations, the priors will have a far greater influence on the resulting estimates. Because there are no priors specified in your Stan model, they default to being uniform priors across the entire range of support.

In other words, your priors are currently:

\mu \sim U(-\infty,\infty)
\sigma \sim U(0,\infty)

This combination of infinitely large priors with a small sample size is likely biasing your estimates.

2 Likes
6

Note that the dataset page explicitly mentions the prior - they have the improper uniform prior on \mu and the impropmer 1/\sigma on \sigma, so you would IMHO need to add something like target += 1/sigma to your model.

1 Like