Why am I getting infinite R-Hat values?

Modeling
1

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?

2

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