Error: Log probability evaluates to log(0), i.e. negative infinity. (Multivariate Gaussian)

I would appreciate it if anyone can help me with this issue, I have changed my priors, however, it does not work. I try to simulate from a multivariate normal distribution based on a given dataset.

The issue is:

Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity
.
.
.
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]]) : Initialization failed.

and my model is:

model<- '
data {
  int<lower=0> N;
  vector[15] x[N];
}

transformed data {
  vector[15] mus = rep_vector(0, 15);
}

parameters {
  cholesky_factor_corr[15] chol;
  vector<lower=0>[15] sigma;
}

transformed parameters {
  matrix[15, 15] chol_cov = diag_pre_multiply(sigma, chol);
}

model {
  chol ~ lkj_corr_cholesky(1);
  sigma ~ cauchy(0, 5);
  x ~ multi_normal_cholesky(mus, chol_cov);
}

generated quantities {
vector[15] x_sim[N];
for(i in 1:N) {
      x_sim[i] = multi_normal_cholesky_rng(mus, chol_cov);
  }
}
'

Okay, I simulated some data from a random covariance matrix in python and used your model to recover it. (NOTE: I used a zero mean for the random multivariate normal draws).

I initialized the above model with inits=0 and it ran fine. Moreover, with N=10000 I recovered my covariance matrix.

Can you explain more what data you are using (i.e. are you sure that your data has zero mean?)

You could augment your model by putting mus in the parameters block and giving it a tight prior around 0. Then if your data doesn’t have zero mean, the model will still have a chance to fit and you would see what the model wants the mean to be.

That is, change your model to:

model<- '
data {
  int<lower=0> N;
  vector[15] x[N];
}

parameters {
  vector[15] mus;
  cholesky_factor_corr[15] chol;
  vector<lower=0>[15] sigma;
}

transformed parameters {
  matrix[15, 15] chol_cov = diag_pre_multiply(sigma, chol);
}

model {
  mus ~ normal(0, 0.1); # Or even 0.01 if you are *really* sure
  chol ~ lkj_corr_cholesky(1);
  sigma ~ cauchy(0, 5);
  x ~ multi_normal_cholesky(mus, chol_cov);
}

[/quote]

Thank you so much for your reply! I also standardized my data, but I got the same error.
There are 50 variables and have small variances, and covariances are so close to zero as well.
Do you think the problem is the covariance chol ~ lkj_corr_cholesky(1)?

model<- '
data {
  int<lower=0> N;
  vector[50] x[N];
}

parameters {
  vector[50] mus;
  cholesky_factor_corr[50] chol;
  vector<lower=0>[50] sigma;
}

transformed parameters {
  matrix[50, 50] chol_cov = diag_pre_multiply(sigma, chol);
}

model {
  mus ~ normal(0, 0.01); 
  chol ~ lkj_corr_cholesky(1);
  sigma ~ cauchy(0, 5);
  x ~ multi_normal_cholesky(mus, chol_cov);
}

generated quantities {
vector[50] x_sim[N];
for(i in 1:N) {
      x_sim[i] = multi_normal_cholesky_rng(mus, chol_cov);
  }
}
'

I highly appreciate your time!

Should I attach my data?

I noticed the model works for the same datasets but by decreasing the number of variables to 20.

With that information, you could try a prior of

chol ~ lkj_corr_cholesky(10)

which would indicate a higher belief that the correlation is the identity.

Still, it is interesting that you are still getting log(0) did you try initializing with all zeros? (I don’t know how to do this in R, but in cmdstanpy you add inits=0 to the sampler)

Thanks so much! will do that and let you know the result. Should I change the number of iterations as the dataset has 50 variables?

The number of iterations is fine. Although for testing, you can reduce them down to a few to make sure things are compiling and running

Good luck!

Thanks! But unfortunately, got the same result when I tried to sample the whole dataset.


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:

I might split the datasets into 3 batches of data.

Try the onion method for that corr matrix Using the onion method in occupancy model hierarchical model - #2 by spinkney

Thanks so much for your reply! I will apply it and let you know the result.