I am getting started with Stan and I’m trying to fit a simple regression model - multivariate normal prior, normal likelihood. However I get an error message saying “The largest R-hat is Inf, indicating chains have not mixed.”

I’ve narrowed it down to the covariance matrix I use in my normal prior. If I take a random covariance matrix the R-Hat values are finite.

The problem matrix is formed in R as:

```
Inv <- function(R) {
chol2inv(chol(R))
}
# Data is in X, this is the formula for the covariance
# matrix for the coefficients fit with ridge regression
xpxi <- Inv(t(X) %*% X + lambda * diag(ncol(X)))
cov_betas <- xpxi %*% t(X) %*% X %*% t(xpxi)
cov_betas_L <- t(chol(cov_betas))
kappa(cov_betas_L) #336649.1, ill-conditioned but not too bad
stan_data <- list(N = nrow(X), M = ncol(X), X = X, y = y,
prior_betas = betas, prior_sigma = cov_betas_L)
```

Stan model:

```
data {
int<lower = 1> N;
int<lower = 1> M;
vector[N] y;
matrix[N, M] X;
matrix[M, M] prior_sigma;
vector[M] prior_betas;
}
parameters {
vector[M] beta;
real<lower = 0> y_sigma;
}
model {
beta ~ multi_normal_cholesky(prior_betas, prior_sigma);
for (n in 1:N) {
target += normal_lpdf(y[n] | beta, y_sigma);
}
}
```

If instead I use the following R code, the R-Hat values are finite.

```
cov_betas_L <- matrix(rnorm(nrow(X) * ncol(X)), nrow = nrow(X))
cov_betas_L <- t(cov_betas_L) %*% cov_betas_L
cov_betas_L <- t(chol(cov_betas_L))
stan_data <- list(N = nrow(X), M = ncol(X), X = X, y = y,
prior_betas = betas, prior_sigma = cov_betas_L)
```

What is wrong with my matrix?