Different results after translation from BUGS to Stan

Modeling
1

I have translated the codes from BUGS, Meyer and Millar (1999) about Bayesian schaefer model, to the Stan.

However, the results are partially different from the original results.

Especially, estimate of K was too big and the shape of a posterior plot had a long tails comparing to the original results of Meyer and Millar (1999)

Could I know the reasons of differences whether from my incorrect translation of coding or difference in algorithm of MCMC between Gibbs sampling and NUTS.

Data file : tuna.dat

C    I
15.9 61.89
25.7 78.98
28.5 55.59
23.7 44.61
25.0 56.89
33.3 38.27
28.2 33.84
19.7 36.13
17.5 41.95
19.3 36.63
21.6 36.33
23.1 38.82
22.5 34.32
22.5 37.64
23.6 34.01
29.1 32.16
14.4 26.88
13.2 36.61
28.4 30.07
34.6 30.75
37.5 23.36
25.9 22.36
25.3 21.91

Stan code file : surplus.stan

data {
  int<lower=0> N; // Number of row
  real<lower=0> C[N]; // Catch
  real<lower=0> I[N]; // CPUE
}
parameters {
  real<lower=0.01, upper=2.0> P[N];
  real<lower=0.01, upper=1.2> r;
  real<lower=10, upper=1000> K;
  real<lower=0.01, upper=0.5> q;
  real sigma2;
  real tau2;
}
transformed parameters {
real Pmed[N];
real Imed[N];
real sigma;
real tau;
real MSY;
real BMSY;
real EMSY;
real P1990;
real B1990;

sigma = sqrt(sigma2);
tau = sqrt(tau2);
MSY = r*K/4;
BMSY = K/2;
EMSY = r/(2*q);
P1990 = P[N]+r*P[N]*(1-P[N]) -C[N]/K;
B1990 = P1990*K;
Pmed[1] = 0;
for (t in 2:N) Pmed[t] = log(P[t-1] + r*P[t-1]*(1-P[t-1]) - C[t-1]/K);
for (t in 1:N) Imed[t] = log(q*K*P[t]);
}

model {
//Priors
K ~ lognormal(5.042905,3.7603664);
r ~ lognormal(-1.38,3.845);
q ~ inv_gamma(0.001,0.001);
sigma2 ~ inv_gamma(3.785518,0.010223);
tau2 ~ inv_gamma(1.708603,0.008613854);

//Joint likelihoods
P[1] ~ lognormal(Pmed[1],sigma);
for (t in 2:N) P[t] ~ lognormal(Pmed[t],sigma);
for (t in 1:N) I[t] ~ lognormal(Imed[t],tau);
}

R source file : surplus.R

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
setwd("C:/Users/KANGHJ/Google 드라이브/Classes/Stan/Practice/MMB")

# Set inputs
df <- read.table("tuna.dat",header = TRUE)
inp <- list(N = length(df$C),
            C = df$C,
            I = df$I)

# Set initial values
ini <- list(P=c(0.99,0.98,0.96,0.94,0.92,0.90,0.88,0.86,0.84,0.82,
                0.80,0.78,0.76,0.74,0.72,0.70,0.68,0.66,0.64,0.62,
                0.60,0.58,0.56),
            r=0.8,
            K=200,
            iq=5,
            sigma2=0.01,
            tau2=0.01)

# Visualize results
set.seed(100)
fit <- stan(file = 'surplus.stan', data = inp, iter = 10000, warmup = 1000, chains= 1, control = list(adapt_delta = 0.99, max_treedepth = 15))
print(fit)


### return an array of three dimensions: iterations, chains, parameters 
a <- extract(fit, inc_warmup = TRUE, permuted = FALSE) 

#bayesplot
library("bayesplot")
library("rstanarm")
library("ggplot2")

plot_title <- ggtitle("Posterior distributions", "with medians and 80% intervals")

color_scheme_set("blue")
mcmc_areas(a, regex_pars = c("K","MSY","BMSY","lp__"), prob = 0.95) + plot_title

color_scheme_set("mix-blue-red")
mcmc_trace(a, pars = c("K","r"), n_warmup = 1000, 
           facet_args = list(nrow = 2, labeller = label_parsed)) + facet_text(size = 15);

My stan code : https://drive.google.com/file/d/0BwKfdyfv6WyPcE5iYkRDSjVNNWc/view?usp=sharing
My data code : https://drive.google.com/file/d/0BwKfdyfv6WyPaWtTN1ZURkIweUU/view?usp=sharing
My rstan code : https://drive.google.com/file/d/0BwKfdyfv6WyPdEFPTkNqUmhoeGc/view?usp=shari

2

That’s a very precise code in some detail:

sigma2 ~ inv_gamma(3.785518,0.010223);

Not so precise, when looking at the lower and upper limits. Also in BUGS codes I see some
truncation, whereas in Stan I don’t see.

3

Translating complex models can always be challenging, but even with a perfect translation you may not get consistent results because the algorithms in BUGS are known to more fragile and prone to bias than those in Stan. We are working are methods to validate BUGS but in general it’s much harder than validating models written in Stan for which we have much more diagnostic information.

4

My version of this statement would be “Stan gives you warnings about divergences and maxtreedepth when it might be possibly biased, BUGS is not able to give such an indication.”. But what he said is correct.

(Ie Stan tells you when it’s broken, BUGS doesn’t, so in the even that Stan isn’t giving any messages and Rhat and ESS are ok, I’d trust Stan.")

5

We often see this and almost always, Stan is providing more accurate results once the same model is being implemented.

The way to test is to generate fake data from the prior and fit it.