Error message: Random variable[1] is nan, but must not be nan!

Modeling
1

Hello everyone,

I am trying to use stan for the first time. The user defined probability density function in my model is

Pr(\boldsymbol{\alpha}) = \frac{\sum_{b'} \exp \left\{ \frac{1} {\delta} \boldsymbol{\alpha}^{T} \mathbf{z}_{b'} + \frac{1}{\delta} \boldsymbol{\omega}^T \mathbf{x}_{b'} - \frac{1}{\delta} P_{b'} \right \}} {\sum_{b'} \exp \left\{ \frac{1}{\delta} \boldsymbol{\omega}^T \mathbf{x}_{b'} - \frac{1}{\delta} P_{b'} + \frac{1}{\delta} \mathbf{z}_{b'}^{T} \boldsymbol{\mu} + \frac{1}{2\delta^2} \mathbf{z}_{b'}^{T} \boldsymbol{\Sigma} \mathbf{z}_{b'} \right\}} f(\boldsymbol{\alpha})
where f(\boldsymbol{\alpha}) is multivariate normal distribution.

And here is my code for the log of the above probability density function:

functions {
  real smvn_lpdf(vector alphas, vector mus, matrix covmats, vector omega, real delta, matrix all_Z, matrix all_X, vector all_P) { 
    
    // Define additional local variables
    vector[256] IMP; 
    real IV; // log of the numerator
    real logK;  // log of the denominator
    real logW; 
    real lprob; 
    
    // determine IMP
    IMP = 1/(delta)*all_X*omega - 1/delta*all_P;

    // determine IV 
    IV = log_sum_exp(1/delta*all_Z*alphas + IMP);

    // calculate logK  
    logK = log_sum_exp(IMP + 1/delta*all_Z*mus + 1/(2*delta^2)*diagonal(all_Z*covmats*all_Z'));

    // calculate logW
    logW = IV - logK;
    
    // calculate the log likelihood 
    lprob = logW + multi_normal_lpdf(alphas | mus, covmats);
    
    // return the log likelihood 
    return lprob;
  }
}

When I ran the model, I got the following error message

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: multi_normal_lpdf: Random variable[1] is nan, but must not be nan!
Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : "
[2] "  c++ exception (unknown reason)"

Could anyone suggest why the random variable is nan, and how I can fix it?

The rest of my code is below:

data {
  int<lower=1> J; 
  int<lower=1> R; 
  int<lower=1> S[R]; 
  int<lower=3> C; 
  int<lower=1> W; 
  int<lower=1> N; 
  int Y[N];  
  matrix[N*3, J-1] Z;  
  matrix[N*3, W] X2; 
  vector[N*3] P; 
  matrix[256, J-1] all_Z;  
  matrix[256, W] all_X;  
  vector[256] all_P; 
}

parameters {
  vector[J-1] mus;
  vector<lower=0>[J-1] sigmas;
  corr_matrix[J-1] Thetas;
  vector[W-1] omegas;
  real<lower = 0> delta;
}

transformed parameters {
  vector[W] omega = append_row(omegas, -sum(omegas));
  cov_matrix[J-1] covmats = quad_form_diag(Thetas, sigmas); 
}

model {
  int yr;
  int xr;
  yr = 1;
  xr = 1;
  
  // hyperpriors
  mus ~ normal(0, 1000);
  sigmas ~ cauchy(0, 2.5); 
  
  Thetas ~ lkj_corr(2);
  
  // priors
  omegas ~ normal(0, 1000);
  delta ~ gamma(0.001, 1000);
  
  // likelihood
  for (r in 1:R) { // for each individual
    int ypos;
    int xpos;
    vector[J-1] alphas;

    ypos = yr;
    xpos = xr;
    alphas ~ smvn(mus, covmats, omega, delta, all_Z, all_X, all_P);	
    
    for (s in 1:S[r]) {  // for each scenario
    // Calculate the choice probability
    Y[ypos] ~ categorical_logit(block(Z, xpos, 1, 2, 6)*alphas
              +block(X, xpos, 1, 2, 27)*omega-1/delta*segment(P, xpos, C)); 
    ypos = ypos + 1;
    xpos = xpos + C;
    }
    yr = yr + S[r];
    xr = xr + S[r]*C;
  }
}

Really appreciate your help!

2

I am a novice in Stan myself, so this might not be helpful. Looking at this part of your code

vector[J-1] alphas; 
ypos = yr; 
xpos = xr; 
alphas ~ smvn(mus, covmats, omega, delta, all_Z, all_X, all_P);

you declare alphas but it is not given any value. You then calculate the log probability of alphas given mus, covmats, omega, delta, all_Z, all_X, all_P. But with alphas not having any value, I would expect that to cause an error. I would expect alphas to be data or transformed data.

1 Like
3

A tilde statement does not assign a value but only computes the probability of a preassigned value. If you want to draw random samples of alpha it needs to be declared in the parameters block.

parameters {
  ...
  real alpha[R,J-1];
}
model {
  ...
  for (r in 1:R) {
    ...
    alpha[r] ~ smvn(..);
    ...
  }
}
3 Likes
4

Thank you for your help! I declared alpha in the parameter block, and it works!

5

Just to be pedantic, alpha[r] ~ smvn(..); isn’t going to draw random samples of alpha in the same way that in R the code:

x = rnorm(1)

draws random numbers. The ~ syntax in Stan increments the log density that the MCMC sampler is working on, which is different (you’ll get posterior samples for alpha given all your data and everything else in the model, which probably won’t have the distribution smvn(..)). Ask if the difference isn’t clear!

Edit: It can be confusing

6

Thank you @bbbales2 for pointing this out. I am not aware of the difference. Maybe I am not careful enough when reading/learning from the online stan examples. Is this difference explained in stan document?

The user defined distribution smvn is the prior distribution for alpha. Thus, isn’t it correct to say that the posterior samples for alpha may not have the smvn distribution?

Thank you for your response!

1 Like
7

Nono, I just brought it up to be pedantic. You’re good.

Yeah, that’s the way to say it. It’s possible they’d have an smvn distribution with different parameters than your priors, but it only works out that way in super simple cases.

The posterior draws from Stan for alpha will be draws from the posterior of the distribution you’re sampling, which is a distribution proportional to your likelihood times prior.

8

Got it! Thanks!