Translating winBUGS code to Stan

Hello all,

I am trying to translate WinBUGS code to Stan, however there seem to be some errors. There may be an obvious mistake I don’t see, since I am new to Stan. Hope that there is someone out there who is kind enough to skim through my code below and enlightened me. Thank you!

Below the WinBUGS code and my Stan version.


model (winBUGS)
{

for (j in 1:30)
{

	# We start by defining that participant j has a unique parameter-value 
	# on each of the six parameters: alpha, beta, gamma, delta, lambda, and luce 
	# (the choice-rule parameter).

	alpha[j] <- phi(alpha.phi[j])
	gamma[j] <- phi(gamma.phi[j])
	luce[j]   <- exp(lluce[j]) 
	
	# We put group-level normal¥s on the individual parameters:

	alpha.phi[j] ~ dnorm(mu.phi.alpha,tau.phi.alpha)I(-3, 3)
	gamma.phi[j] ~ dnorm(mu.phi.gamma,tau.phi.gamma)I(-3, 3)
        	lluce[j]   ~ dnorm(lmu.luce, ltau.luce)
}

# Here priors for the hyperdistributions are defined:
mu.phi.alpha ~ dnorm(0,1)
tau.phi.alpha <- pow(sigma.phi.alpha,-2)
sigma.phi.alpha ~ dunif(0,10)

mu.phi.gamma ~ dnorm(0,1)
tau.phi.gamma <- pow(sigma.phi.gamma,-2)
sigma.phi.gamma ~ dunif(0,10)

  	lmu.luce    ~ dunif(-2.3, 1.61)   
  	ltau.luce  <- pow(lsigma.luce,-2)
  	lsigma.luce ~ dunif(0,1.13)       


# To obtain the mean of the hyper distribution on the wanted scale:
mu.alpha <- phi(mu.phi.alpha) 
mu.gamma <- phi(mu.phi.gamma) 
  	mu.luce   <- exp(lmu.luce)


for (j in 1:30)# Subject-loop
{
	for (i in 1:60)# Item-Loop, positive gambles,gamble A
	{
 			v.x.a[i,j] <- pow(prospects.a[i,1],alpha[j])  
  			w.x.a[i,j] <- pow(prospects.a[i,4],gamma[j]) / pow(z.a[i,j],(1/gamma[j])) 
  			z.a[i,j]   <- pow(prospects.a[i,2],gamma[j]) + pow(prospects.a[i,4],gamma[j]) 
  			Vf.a[i,j]  <- w.x.a[i,j] * v.x.a[i,j]
	}


	for (i in 1:60)# Item-Loop, positive gambles,gamble B
	{
  			v.x.b[i,j] <- pow(prospects.b[i,1],alpha[j])        	 
  			w.y.b[i,j] <- pow(prospects.b[i,4],gamma[j]) / pow(z.b[i,j],(1/gamma[j])) 
  			z.b[i,j]   <- pow(prospects.b[i,2],gamma[j]) + pow(prospects.b[i,4],gamma[j]) 
  			Vf.b[i,j]  <- w.x.b[i,j] * v.x.b[i,j]
  
	}

	for (i in 1:60)# Item-Loop, positive gambles,choice-rule
	{
  
		binval[i,j] <- (1)/(1+exp((-1*luce[j])*(Vf.b[i,j]-Vf.a[i,j])))
		rawdata[i,j] ~ dbern(binval[i,j])
	}

// Model: cumulative prospect theory.
// Takes choice data and gamble attributes as input
data {
  int<lower=0> N;
  int<lower=0> Nsubs;
  real prospects[31, 4]; // table gamble specifications
  int <lower=0,upper=1> choicedata[N, Nsubs]; // choice outcome, 0 or 1
}

parameters {
  
  real mu_phi_alpha;
  real mu_phi_gamma;
  real lmu_luce ;
  
  real alpha_phi[Nsubs];
  real gamma_phi[Nsubs];
  real lluce[Nsubs];
  
  real <lower = 0, upper = 10> sigma_phi_alpha;
  real <lower = 0, upper = 10> sigma_phi_gamma;
  real <lower = 0, upper = 1.13> lsigma_luce;
}

transformed parameters {
  
  real <lower = 0, upper = 2>  alpha[Nsubs]; 
  real <lower = 0, upper = 1>  gamma[Nsubs];
  real <lower = 0, upper = 10>  luce[Nsubs];

  real mu_alpha;
  real mu_gamma;
  real mu_luce ;
  
  alpha = Phi(alpha_phi);
  gamma = Phi(gamma_phi);
  luce = exp(lluce);

  mu_alpha = Phi(mu_phi_alpha);
  mu_gamma = Phi(mu_phi_gamma);
  mu_luce = exp(lmu_luce);
}

model { 
  real vxa[N, Nsubs];
  real wxa[N, Nsubs];
  real vxb[N, Nsubs];
  real wxb[N, Nsubs];
  real za[N, Nsubs];
  real zb[N, Nsubs];
  real Vfa[N, Nsubs];
  real Vfb[N, Nsubs];
  real diff[N, Nsubs];

  // target += normal_lpdf(mu_phi_alpha|0,1);
  // target += normal_lpdf(mu_phi_gamma|0,1);
  // target += uniform_lpdf(lmu_luce|-2.3,1.61);
  
  mu_phi_alpha ~ normal(0, 1);
  mu_phi_gamma ~ normal(0, 1);
  lmu_luce ~ uniform(-2.3, 1.61);

  // target += uniform_lpdf( sigma_phi_alpha|0,10);
  // target += uniform_lpdf( sigma_phi_gamma|0,10);
  // target += uniform_lpdf( lsigma_luce|0,1.13);
  
  sigma_phi_alpha ~ uniform(0, 10);
  sigma_phi_gamma ~ uniform(0, 10);
  lsigma_luce ~ uniform(0, 1.13);

  for (j in 1:Nsubs) { 
    
    // target += normal_lpdf(alpha_phi[j]|mu_phi_alpha, sigma_phi_alpha);
    // target += normal_lpdf(gamma_phi[j]|mu_phi_gamma, sigma_phi_gamma);
    // target +=  normal_lpdf(lluce[j]|lmu_luce, lsigma_luce);

    alpha_phi[j]~ normal(mu_phi_alpha, sigma_phi_alpha);
    gamma_phi[j] ~ normal(mu_phi_gamma, sigma_phi_gamma);
    lluce[j] ~ normal(lmu_luce, lsigma_luce);
  } 

  for (j in 1:Nsubs)// Subject-loop
  {
    for (i in 1:N)// Item-Loop, positive gambles,gamble A
    {
      vxa[i,j] = pow(prospects[i,1], alpha[j]);
      za[i,j]  = pow(prospects[i,2],gamma[j]) + pow(1-prospects[i,2],gamma[j]);
      wxa[i,j] = pow(prospects[i,2],gamma[j]) / pow(za[i,j],(1/gamma[j])) ; 
      Vfa[i,j]  =wxa[i,j] * vxa[i,j];
      // gamble B
      vxb[i,j] = pow(prospects[i,3],alpha[j])   ;  
      zb[i,j]  = pow(prospects[i,4],gamma[j]) + pow(1-prospects[i,4],gamma[j]) ;
      wxb[i,j] = pow(prospects[i,4],gamma[j]) / pow(zb[i,j],(1/gamma[j])) ;
      Vfb[i,j] = wxb[i,j] * vxb[i,j];
      diff[i,j] = (Vfb[i,j]-Vfa[i,j]);
    }
    
    for (i in 1:N) //Item-Loop, positive gambles,choice-rule
     {
    choicedata[i,j] ~ bernoulli_logit((1)/(1+exp((-1*luce[j])*(Vfb[i,j]-Vfa[i,j]))));
     }
  }
}




Paper describing code: https://www.ejwagenmakers.com/2011/NilssonEtAl2011.pdf

(I keep getting following error message:

  1. Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable. If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform. Left-hand-side of sampling statement: luce[j] ~ normal(…)
    )

The warning (not an error) says that you have a line that attempts to place a prior on luce, which would require a Jacobian adjustment. However, in the code that you show here, you don’t place a prior on luce but rather on lluce. This should not require a Jacobian adjustment. I wondered whether this resulted from clumsy parsing, but it seems not to. The following program compiles and samples without warning in both rstan 2.21 and 2.26

parameters {
  real lluce;
}
transformed parameters {
  real luce = lluce^2;
}
model { 
   lluce ~ std_normal();
}

But the following program (correctly) gives this Jacobian warning in both versions of rstan

parameters {
  real lluce;
}
transformed parameters {
  real luce = lluce^2;
}
model { 
   luce ~ std_normal();
}

Is it perhaps the case that you received this warning from a slightly different Stan program than the one you included in your post? Alternatively, do you still get the warning if you run the first Stan program above?

1 Like

Thank you for your response! I ran the program again, this time the warning did not appeared. It must have been the warning from a slightly different Stan program as you suggested. Thanks for catching that.

Now initialisation fails when I run the program with simulated data. (This has happened previously, but I thought it was a result of the previous warning message regarding the Jacobian adjustment. See Exception: bernoulli_logit_lpmf: Logit transformed probability parameter is nan, but must not be nan!):

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: bernoulli_logit_lpmf: Logit transformed probability parameter is nan, but must be not nan! (in ‘model4ab459f32f4_cpt_2’ at line 101)

Any ideas what the problem might be? Thanks again!

Don’t write 1/(1+exp(-1*...)) when Stan has inv_logit(...)

    choicedata[i,j] ~ bernoulli_logit(inv_logit(luce[j]*(Vfb[i,j]-Vfa[i,j])));

And actually the _logit in bernoulli_logit already takes care of that

    choicedata[i,j] ~ bernoulli_logit(Vfb[i,j]-Vfa[i,j]);

This is common with BUGS but Stan does not require declaring all the variables at the top level. You can just declare the variable inside the loop where it is used, for instance

  // no toplevel arrays needed
  for (j in 1:Nsubs)// Subject-loop
  {
    for (i in 1:N)// Item-Loop, positive gambles,gamble A
    {
      real vxa = pow(prospects[i,1], alpha[j]);
      real za  = pow(prospects[i,2],gamma[j]) + pow(1-prospects[i,2],gamma[j]);
      real wxa = pow(prospects[i,2],gamma[j]) / pow(za,(1/gamma[j])) ; 
      real Vfa = wxa * vxa;
      // gamble B
      real vxb = pow(prospects[i,3],alpha[j])   ;  
      real zb  = pow(prospects[i,4],gamma[j]) + pow(1-prospects[i,4],gamma[j]) ;
      real wxb = pow(prospects[i,4],gamma[j]) / pow(zb,(1/gamma[j])) ;
      real Vfb = wxb * vxb;
      real diff = (Vfb-Vfa);
      // positive gambles,choice-rule
      choicedata[i,j] ~ bernoulli_logit(luce[j]*diff);
    }
  }

Parameters with uniform distribution should be declared with explicit bounds (as you have correctly done with sigma_phi_alpha, sigma_phi_gamma and lsigma_luce)

parameters {
  real<lower=-2.3, upper=1.61> lmu_luce ;

That suggest something goes wrong when computing the probability. The most obviously suspect…are you sure prospects[i,n] is always between 0 and 1? Your BUGS model does not have a pow(1-prospects[i,2], gamma) term in it.

1 Like

Thank you, this is brilliant advice! I adapted my code accordingly, unfortunately my RStudio keeps crashing when running the model. It also crashes when I run simple models, or the examples provided by Stan. Any idea why that happens? Thanks again!

One reason that Rstudio crashes these days is if you have mismatched versions of packages rstan and StanHeaders. There’s no guarantee that this is the problem, but if it is you should be able to fix by running (in a fresh R session) one of the following:

To use a more up-to-date version of rstan

remove.packages(c("StanHeaders", "rstan"))
install.packages("StanHeaders", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

To use the CRAN version

remove.packages(c("StanHeaders", "rstan"))
install.packages("rstan")
1 Like

Thank you, works perfectly now!