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:
-
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(…)
)