Error message with Betancourt specification of horseshoe prior

Modeling
1

Hi all,

I’m trying to run an analysis using the horseshoe specification in Bayes Sparse Regression (2018). The full code is here

library(rstan)
library(loo)
library(bayesplot)
library(xtable)

## Read in data ##
PISA18sampleScale <- read.csv(file.choose(),header=T)

n=nrow(PISA18sampleScale)
X = PISA18sampleScale[,2:11]
readscore = PISA18sampleScale[,1]
p=ncol(X)

data.list <- list(n=n, p=p, X=X, readscore=readscore)

modelString = "
data {
  int<lower=1> n; // Number of data
  int<lower=1> p; // Number of covariates
  matrix[n, p] X;
  real readscore[n];
}

parameters {
  vector[p] beta_tilde;
  vector<lower=0>[p] lambda;
  real<lower=0> tau_tilde;
  real alpha;
  real<lower=0> sigma;
}

model {
  // tau ~ cauchy(0, sigma)
  // beta ~ normal(0, tau * lambda)
  vector[p] beta = beta_tilde .* lambda * sigma * tau_tilde;
  
  beta_tilde ~ normal(0, 1);
  lambda ~ cauchy(0, 1);
  tau_tilde ~ cauchy(0, 1);
  
  alpha ~ normal(0, 1);
  sigma ~ cauchy(0, 1);
  
  readscore ~ normal(X' * beta + alpha, sigma);
}
"
# Start estimation

nChains = 4
nIter= 10000
thinSteps = 10
warmupSteps = floor(nIter/2)

readscore = data.list$readscore

myfitHorseshoe = stan(data=data.list,model_code=modelString,chains=nChains,
                  iter=nIter,warmup=warmupSteps,thin=thinSteps)

I’m getting the following error message

SAMPLING FOR MODEL '6a7711900e92e051ef7d7e1268a95faf' NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: multiply: Columns of m1 (100) and Rows of m2 (10) must match in size  (in 'model98d60f2608e_6a7711900e92e051ef7d7e1268a95faf' at line 29)

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                      
[2] "  Exception: multiply: Columns of m1 (100) and Rows of m2 (10) must match in size  (in 'model98d60f2608e_6a7711900e92e051ef7d7e1268a95faf' at line 29)"
error occurred during calling the sampler; sampling not done
>

In addition this error, I’m curious why the specification tau and beta are commented out.

2

I take it your code is this with M renamed to p. You have matrix X transposed, Betancourt defines matrix[M,N] X but your code says matrix[n,p] X.
You can leave that and simply remove the transpose operator ' from the last line

  readscore ~ normal(X * beta + alpha, sigma); // not X'

The model is parametrized with tau_tilde and beta_tilde instead of (more natural?) tau and beta. (I guess it’s like noncentered parametrization or something?) The model block defines the distributions for these parameter and the comments document the implied distributions for the derived parameters.
In other words, tau = tau_tilde * sigma so then tau_tilde ~ cauchy(0, 1) is an indirect way to declare tau ~ cauchy(0, sigma).

3

Thanks very much. That was a silly mistake. It seems to be working but throwing divergent transition errors. I’m playing with the adapt_delta value and hopefully that will help.