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.