Transformed parameter evaluates to #QNAN

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);

}

Probably overflows. You should be using bernoulli_logit_lpmf in line 41 to avoid this problem. Also, you should not be double looping. Do

for (n in 1:N) {
  target += benoulli_logit_lpmf(rawdata[n,] | ...);
}

but this may require that you use vectors for some things rather than multidimensional real arrays.

Thanks! After making the recommended changes, I still get the same error, unfortunately. I was able to figure out that it has something to do with the code that multiplies 0 by -1. When I remove observations that would result in -1 * 0, the the NAN error goes away (although I have another error in its place that I’m working to resolve). Has anyone seen this issue before?

Have you tried printing the intermediate computations to make sure it’s what you think it is? It’s a pretty complicated transformed parameters block.

A quick way to check is to comment out your model block and throw normal(0, 1) priors on all your parameters. That should sample just fine and you can focus on figuring out what your transformed parameters block is actually doing.

And… no. -1 * 0 shouldn’t result in NaN.

Thanks for the suggestions! I figured out that the issue has something to do with this part of the code:

pow((-1*prospecta[n,1]), beta[j])

where beta is constrained to be in [0,1] and prospect[n,1] is in (-inf, 0]. For whatever reason, this function returns NAN. When I change the data set so that prospect[n,1] is positive (so there’s no need to multiply by -1), the code runs just fine.

This is the n-th root of a negative number, no? If so, well, that’s a pretty good reason to get a NaN, imo.

Just for future reference when troubleshooting, you can do something like:

if (!is.finite(pow((-1*prospecta[n,1]), beta[j]))) {
  print(prospecta[n,1]):
  print(beta[j]);
}

Then you’re not left with “for whatever reason” as the answer.

I guess I didn’t explain it so well. In the data set, prospect[n,1] is negative. The formula is supposed to flip the sign before taking the n-th root.

OK, but what do you mean exactly by “flipping the sign”? If prospect[n,1] is negative, you can always do abs(prospect[n,1]). But if prospect[n,1] > 0, flipping its sign will run into the same problem. Are you sure the variable is always negative? @sakrejda’s advice is also sound, imho.

@maxbiostat, thanks a lot for your help. I used @sakrejda’s advice so I could pinpoint the problem.

print(n)
print(prospect[n,3]);
print(-1*prospect[n,3]);
print(pow((-1*prospecta[n,3]), beta[j])); 

reveals:

67
0
-0
-1.#IND

So Stan appears to think that -0 is a negative number, and outputs an error because the n-th root of a negative number is not real.

2 Likes

Not Stan so much as IEEE arithmetic. We’re just delegating these operations to the standard C++ libraries.

Now that’s nasty. I don’t know how to fix this in general in our library without an expensive zero normalization step everywhere. Is there a function of ours that’s checking for negativity and failing? That we could probalby control.