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

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

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

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


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

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

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

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

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?