Hi all,

I’m having some trouble with this hierarchical non-centered linear model I’m attempting to run. Whenever I begin sampling, I get the error: Log probability evaluates to log(0), i.e. negative infinity.

Stan can’t start sampling from this initial value.

I’ve run this model on simulated data with the same result, and I friend of mine attempting to use the same model is also receiving this error. This leads us to believe that there’s potentially some sort of misspecification in the model, but we’re not sure where or what the issue would be. Any help in the right direction would be super helpful.

In terms of data, both my friend and I have outcome variables that are constrained to be positive: mine is an index between 0 and 100 (though I’ve tried transforming the data on an offset log scale before running the model with the same result), and my friend is analyzing sales data with the outcome between $0 and $2million.

Code:

```
// Index values, observations, and covariates.
data {
int<lower = 1> N; // Number of observations.
int<lower = 1> K; // Number of groups/clusters.
int<lower = 1> I; // Number of observation-level covariates.
vector[N] y; // Vector of observations.
array[N] int<lower = 1, upper = K> g; // Vector of group assignments.
matrix[N, I] X; // Matrix of observation-level covariates.
}
// Parameters.
parameters {
vector[I] gamma; // Vector of population-level parameters.
cholesky_factor_corr[I] L_Omega; // Population model correlation matrix Cholesky factor.
vector[I] tau; // Population model scale parameters.
matrix[I, K] Delta; // Matrix of non-centered observation-level parameters.
real<lower = 0> sigma; // Variance of the observation model.
}
// Transformed parameters.
transformed parameters {
matrix[K, I] V; // Matrix of rescaled non-centered observation-level parameters.
matrix[K, I] Beta; // Matrix of observation-level parameters.
// Compute the non-centered parameterization.
V = (diag_pre_multiply(tau, L_Omega) * Delta)';
for (i in 1:I) {
Beta[,i] = gamma[i] + V[,i];
}
}
// Regression model.
model {
// Hyperpriors.
gamma ~ normal(0, 5);
L_Omega ~ lkj_corr_cholesky(4);
tau ~ normal(0, 5);
// Prior.
sigma ~ normal(0, 5);
// Population model.
to_vector(Delta) ~ normal(0, 1);
// Observation model.
for (n in 1:N) {
y[n] ~ lognormal(X[n,] * Beta[g[n],]', sigma);
}
}
// Generate quantities at each draw.
generated quantities {
// Compute the log-likelihood.
vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = lognormal_lpdf(y[n] | X[n,] * Beta[g[n],]', sigma);
}
// Compute the correlation matrix Omega.
corr_matrix[I] Omega;
Omega = multiply_lower_tri_self_transpose(L_Omega);
}
```

Thanks again for any help you can provide! We’re pretty stumped on this one.