Assistance to convert OpenBugs Code into Stan Code

Modeling
1

Hi all,

I was introduced to OpenBUGS by a colleague and has been using it since. Recently I came across STAN and thought maybe I could give it a try. However there seems to be some error. Is there anyone who is kind enough to vet through my code below and enlightened me. Many thanks.

model{
mu ~ dnorm(0,1.0E-6)
beta ~ dgamma(1.0E-5,1.0E-5)

for (i in 1:6){
delta[i] ~ dnorm(mu,beta)

prec[i] <- 1/(u[i] * u[i])
x[i] ~ dnorm(delta[i],prec[i])

pred[i] ~ dnorm(mu,prec[i])
DOE[i] <- x[i] - pred[i]
}
}

list(
x = c(34.30,32.90,34.53,32.42,31.90,35.80),
u = c(1.03,0.69,0.83,0.29,0.40,0.38)
)

list(beta=1)

I am trying to convert the above into:

data {
  int<lower=0> N;         
  real x[N];            
  real<lower=0> u[N];  
}

parameters {
  real mu;
  real <lower=0> beta;
  real <lower=0> delta[N];
  real <lower=0> pred[N];
}

transformed parameters{
  real uu[N];
  uu[N] = u[N] * u[N];
}

model {
  mu ~ normal(0,1E-6);  //prior for mu
  beta ~ gamma(1E-5,1E-5);  //prior for beta
  delta[N]~normal(mu,beta);
  x[N] ~ normal(delta[N],uu[N]);
  pred[N]~normal(mu,uu[N]);
}

generated quantities{
  real DOE[N];
  DOE[N] = x[N] -pred[N];
}

I saved the above as STAN extension and calling it in R with the following:
```R script
library(rstan)
library(shinystan)

options(mc.cores = 4)
data <-list(N=6, y = c(34.3,32.9,34.53,32.42,31.9,35.8), u=c(1.03,0.69,0.83,0.29,0.4,0.38))

model1 <- stan(file='Openbugs convert HB REM1.stan', data=data)
print(model1)
2

Hi Benny, welcome to Stan!

There are a few things here. Firstly, I’d highly recommend reading the User’s Guide chapter on transitioning from BUGS: https://mc-stan.org/docs/2_25/stan-users-guide/stan-for-bugs-appendix.html

Now, onto your model:

Vectors

When working with multiple values for a single parameter, it is more efficient to use vectors, rather than arrays of real, see this section of the manual for more information on the various data types: https://mc-stan.org/docs/2_25/reference-manual/data-types-chapter.html.

So instead of:

real<lower=0> u[N];

Use:

vector<lower=0>[N] u;

Indexing

Something to be very careful with your models is how you use indexes. Take the following statement for example:

delta[N]~normal(mu,beta);

This is not putting a prior on the whole vector of delta values, only on the sixth value. This is because N is passed as the single integer 6, which means that statement is equivalent to:

delta[6]~normal(mu,beta);

If you want the prior (or any function) to apply to the whole vector, then you can just specify it without any indexing:

delta ~ normal(mu,beta);

Normal Distribution

Just a quick one. Stan parameterises the normal distribution using the mean and standard deviation, not the variance, so you don’t need to square the SD parameter (u)

Predicted Values

To generate the model-predicted values, you should use the _rng function of your desired distribution in the generated quantities block, rather than declaring them as a new parameter.

So rather than having in the model block:

pred[N]~normal(mu,uu[N]);

You would instead have in the generated quantities:

pred = normal_rng(mu, u);

The Fixed Model

I’ll post a corrected version of your model below, but make sure you take some time to see the changes I’ve made and why, and I’d also recommend reading the Stan documentation here: https://mc-stan.org/users/documentation/ to get more familiar with the language.

data {
  int<lower=0> N;         
  vector[N] x;            
  vector<lower=0>[N] u;  
}

parameters {
  real mu;
  real<lower=0> beta;
  vector<lower=0>[N] delta;
}


model {
  mu ~ normal(0,1E-6);  //prior for mu
  beta ~ gamma(1E-5,1E-5);  //prior for beta
  delta ~ normal(mu,beta);
  x ~ normal(delta,u);
}

generated quantities{
  real pred[N];
  vector[N] DOE;

  pred = normal_rng(mu, u);

  // The vectorised rng returns a real[],
  // need to convert to vector for subtraction
  DOE = x - to_vector(pred);
}
5 Likes
3

That also affects the prior parametrization. It should be

  mu ~ normal(0,1000);  //prior for mu
2 Likes
4

Oh good catch, I missed that! Except the original prior was normal(0,1E-6), so the SD prior would be normal(0, 0.001)

5

No, OpenBUGS uses precision instead of variance. You can see this in the original model where it says

prec[i] <- 1/(u[i] * u[i])
x[i] ~ dnorm(delta[i],prec[i])
1 Like
6

Oh right, I was assuming that Benny had already put it in the variance scale, like with the uu parameter, but I should have checked that

1 Like
7

Hi @andrjohns @nhuurre, Great! Many thanks for taking your time to understand my OpenBUGS code and converting to STAN (also at the same time guiding me along!). To be honest, I am a chemist by training, I picked up the stats while on the job and recently tried exploring Bayesian stats. Hope I can learn from you all on this forum.
Thanks again, greatly appreciate your assistance.

3 Likes