# Why am I getting infinite R-Hat values?

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?

Not sure, but maybe change the model to:

``````parameters {
vector[M] beta_z;
real<lower = 0> y_sigma;
}

model {
vector[M] beta = prior_betas + prior_sigma * beta_z;
...
}
``````

That’ll imply the same prior on beta and maybe it’ll be more stable or whatnot.

2 Likes