I’m trying to estimate the parameters of Cumulative Prospect Theory. This requires several transformation be performed to the data and model parameters. When I try run samples, I get the following error:
Exception: bernoulli_lpmf: Probability parameter is 1.#QNAN, but must
be finite! (in ‘C:/Users/mstern10/Documents/HBCPT.stan’ at line 141)
Here’s my code:
data {
int<lower=0> N; //Number of observations
int<lower=0> J; //Number of subjects
real prospecta[N,4]; //Gamble A’s parameters
real prospectb[N,4]; //Gamble B’s parameters
int rawdata[N,J]; //Choice data from subjects
}
parameters {
real<lower=-3, upper=3> alphaphi[J];
real<lower=-3, upper=3> betaphi[J];
real<lower=-3, upper=3> gammaphi[J];
real<lower=-3, upper=3> deltaphi[J];
real lluce[J];
real llambda[J];
real mualphaphi;
real mubetaphi;
real mugammaphi;
real mudeltaphi;
real lmuluce;
real lmulambda;
real<lower=0, upper=10> sigmaalphaphi;
real<lower=0, upper=10> sigmabetaphi;
real<lower=0, upper=10> sigmagammaphi;
real<lower=0, upper=10> sigmadeltaphi;
real<lower=0, upper=1.13> lsigmaluce;
real<lower=0, upper=1.13> lsigmalambda;
}
transformed parameters {
real mualpha;
real mubeta;
real mugamma;
real mudelta;
real muluce;
real mulambda;
real alpha[J];
real beta[J];
real gamma[J];
real delta[J];
real luce[J];
real lambda[J];
real valuea1[N,J];
real valuea2[N,J];
real valueb1[N,J];
real valueb2[N,J];
real za[N,J];
real zb[N,J];
real weighta1[N,J];
real weighta2[N,J];
real weightb1[N,J];
real weightb2[N,J];
real Valuea[N,J];
real Valueb[N,J];
real binval[N,J];
//Transform parameters so they can be used to calculate values;
for (j in 1:J){
alpha[j] = normal_cdf(alphaphi[j], 0, 1);
beta[j] = normal_cdf(betaphi[j], 0, 1);
gamma[j] = normal_cdf(gammaphi[j], 0, 1);
delta[j] = normal_cdf(deltaphi[j], 0, 1);
luce[j] = exp(lluce[j]);
lambda[j] = exp(lambda[j]);
}
for (j in 1:J){
//Calculate the value of positive gambles;
for (n in 1:60){
valuea1[n,j] = pow(prospecta[n,1], alpha[j]);
valuea2[n,j] = pow(prospecta[n,3], alpha[j]);
za[n,j] = pow(prospecta[n,2], gamma[j]) + pow(prospecta[n,4], gamma[j]);
weighta1[n,j] = pow (prospecta[n,2], gamma[j]) / pow(za[n,j],(1/gamma[j]));
weighta2[n,j] = pow (prospecta[n,4], gamma[j]) / pow(za[n,j],(1/gamma[j]));
Valuea[n,j] = weighta1[n,j] * valuea1[n,j] + weighta2[n,j] * valuea2[n,j];
valueb1[n,j] = pow(prospectb[n,1], alpha[j]);
valueb2[n,j] = pow(prospectb[n,3], alpha[j]);
zb[n,j] = pow(prospectb[n,2], gamma[j]) + pow(prospectb[n,4], gamma[j]);
weightb1[n,j] = pow (prospectb[n,2], gamma[j]) / pow(zb[n,j],(1/gamma[j]));
weightb2[n,j] = pow (prospectb[n,4], gamma[j]) / pow(zb[n,j],(1/gamma[j]));
Valueb[n,j] = weightb1[n,j] * valueb1[n,j] + weightb2[n,j] * valueb2[n,j];
}
//Calculate the value of negative gambles (note, original code appears to make error and ignore lambda);
for (n in 61:120){
valuea1[n,j] = (-1) * pow((-1 * prospecta[n,1]), beta[j]);
valuea2[n,j] = (-1) * pow((-1 * prospecta[n,3]), beta[j]);
za[n,j] = pow(prospecta[n,2], delta[j]) + pow(prospecta[n,4], delta[j]);
weighta1[n,j] = pow (prospecta[n,2], delta[j]) / pow(za[n,j],(1/delta[j]));
weighta2[n,j] = pow (prospecta[n,4], delta[j]) / pow(za[n,j],(1/delta[j]));
Valuea[n,j] = weighta1[n,j] * valuea1[n,j] + weighta2[n,j] * valuea2[n,j];
valueb1[n,j] = (-1) * pow((-1 * prospectb[n,1]), beta[j]);
valueb2[n,j] = (-1) * pow((-1 * prospectb[n,3]), beta[j]);
zb[n,j] = pow(prospectb[n,2], delta[j]) + pow(prospectb[n,4], delta[j]);
weightb1[n,j] = pow (prospectb[n,2], delta[j]) / pow(zb[n,j],(1/delta[j]));
weightb2[n,j] = pow (prospectb[n,4], delta[j]) / pow(zb[n,j],(1/delta[j]));
Valueb[n,j] = weightb1[n,j] * valueb1[n,j] + weightb2[n,j] * valueb2[n,j];
}
//Calculate the value of mixed gambles;
for (n in 121:180){
valuea1[n,j] = pow(prospecta[n,1], alpha[j]);
valuea2[n,j] = (-1 * lambda[j]) * pow((-1 * prospecta[n,3]), beta[j]);
za[n,j] = pow(prospecta[n,2], delta[j]) + pow(prospecta[n,4], delta[j]);
weighta1[n,j] = pow (prospecta[n,2], delta[j]) / pow(za[n,j],(1/delta[j]));
weighta2[n,j] = pow (prospecta[n,4], delta[j]) / pow(za[n,j],(1/delta[j]));
Valuea[n,j] = weighta1[n,j] * valuea1[n,j] + weighta2[n,j] * valuea2[n,j];
valueb1[n,j] = pow(prospectb[n,1], alpha[j]);
valueb2[n,j] = (-1 * lambda[j]) * pow((-1 * prospectb[n,3]), beta[j]);
zb[n,j] = pow(prospectb[n,2], gamma[j]) + pow(prospectb[n,4], delta[j]);
weightb1[n,j] = pow (prospectb[n,2], gamma[j]) / pow(zb[n,j],(1/gamma[j]));
weightb2[n,j] = pow (prospectb[n,4], delta[j]) / pow(zb[n,j],(1/delta[j]));
Valueb[n,j] = weightb1[n,j] * valueb1[n,j] + weightb2[n,j] * valueb2[n,j];
}
}
//Calculate the probability of choosing gamble b;
for (j in 1:J){
for (n in 1:N){
binval[n,j] = (1)/(1+exp((-1luce[j])(Valueb[n,j]-Valuea[n,j])));
}
}
mualpha = normal_cdf(mualphaphi, 0, 1);
mubeta = normal_cdf(mubetaphi, 0, 1);
mugamma = normal_cdf(mugammaphi, 0, 1);
mudelta = normal_cdf(mudeltaphi, 0, 1);
mulambda = exp(lmulambda);
muluce = exp(lmuluce);
}
model {
for (j in 1:J) {
for(n in 1:N) {
rawdata[n,j] ~ bernoulli(binval[n,j]);
}
}
for (j in 1:J) {
alphaphi[j] ~ normal(mualphaphi, sigmaalphaphi)T[-3,3];
betaphi[j] ~ normal(mubetaphi, sigmabetaphi)T[-3,3];
gammaphi[j] ~ normal(mugammaphi, sigmagammaphi)T[-3,3];
deltaphi[j] ~ normal(mudeltaphi, sigmadeltaphi)T[-3,3];
lluce[j] ~ normal(lmuluce, lsigmaluce);
llambda[j] ~ normal(lmulambda, lsigmalambda);
}
mualphaphi ~ normal(0,1);
sigmaalphaphi ~ uniform(0,10);
mubetaphi ~ normal(0,1);
sigmabetaphi ~ uniform(0,10);
mugammaphi ~ normal(0,1);
sigmagammaphi ~ uniform(0,10);
mudeltaphi ~ normal(0,1);
sigmadeltaphi ~ uniform(0,10);
lmulambda ~ uniform(-2.3, 1.61);
lsigmalambda ~ uniform(0, 1.13);
lmuluce ~ uniform(-2.3, 1.61);
lsigmaluce ~ uniform(0, 1.13);
}