Hyper prior to obtain sparsity

This is some code about latent dimension variable. The code get stuck for the following error: “Exception thrown at line 69: []: accessing element out of range. index 5 out of range; expecting index to be between 1 and 4; index position = 2Pi”. I know the problem is on the highlighted part below. So what is exactly the problem?

I read the suggestion from avehtari to use horseshoe prior. But before I switch to that, does Stan not able to deal with exponential distribution as a prior?

Simulated data
######## Set up parameters #####
p=10; q=9;D=p+q;
k=4; k1=2; k2=1; K=k+k1+k2;
n=100;
##################### Generate Z ######################
mu_z <- rep(0,K)
Sigma_z <- matrix(0, nrow=K, ncol=K) + diag(K)*1

Z <- mvrnorm(n, mu_z, Sigma_z)
##################### Generate H #####################
H <- matrix(nrow=D, ncol=K)

First k, shared

set.seed(1)############################################ random generator
for(j in 1:k)
{H[,j] = rbinom(D, 1, 2/(j+2))}

k+1 to k+k1, x specific

for(j in (k+1):(k+k1))
{ for(i in 1:p) H[i,j]=1
for(i in (p+1):D) H[i,j]=0 }

k+k1+1 to K, y specific

for(j in (k+k1+1):K)
{ for(i in 1:p) H[i,j] = 0
for(i in (p+1):D) H[i,j] = 1}

Generate E

set.seed(1)############################################ random generator
E_mux1 <- runif(K, min =.01, max = 1)
E_mux <- sort(E_mux1, decreasing = TRUE)
set.seed(1)
e_x <- runif(K, min = .5, max = 5)
E_sigmax <- diag(K)*e_x
E_sigmax
set.seed(1)############################################ random generator
E_x <- mvrnorm(p, E_mux, E_sigmax)
set.seed(1)############################################ random generator
E_muy1 <- runif(K, min =5, max = 40)
E_muy <- sort(E_muy1, decreasing = TRUE)
set.seed(1)############################################ random generator
e_y <- runif(K, min = .5, max = 20)
E_sigmay <- diag(K)*e_y
set.seed(1)############################################ random generator
E_y <- mvrnorm(q, E_muy, E_sigmay)
E <- rbind(E_x, E_y)

Generate Phi

Phi <- E*H

Generate Theta

Theta <- Z%*%t(Phi)
Theta_x <- Theta[,1:10]
Theta_y <- Theta[,11:19]
#######################Generate Data###########################

X

Generate X related parameter

X_p <- matrix(nrow=n, ncol=p,)
{for(j in 1:p)
{ for(i in 1:n)
{X_p[i,j] = exp(Theta_x[i,j])/(1+exp(Theta_x[i,j]))}}}
set.seed(1)############################################ random generator
X<- matrix(nrow=n, ncol=p)
{for(j in 1:p)
for(i in 1:n) {X[i,j] = rbinom(1, 2, X_p[i,j])}}
table(X)

Generate Y

Generate Y related parameter

Y_mu <- matrix(nrow=n, ncol=q)
{for(j in 1:q)
for(i in 1:n) {Y_mu[i,j] = Theta_y[i,j]}}

identical(Y_mu, Theta_y)

set.seed(1)
Y<- matrix(nrow=n, ncol=q)
{for(j in 1:q)
{ for(i in 1:n)
{Y[i,j] = rnorm(1,Theta_y[i,j], sd=1/2)}}}

stan code

library(rstan)
library(parallel)
library(coda)
SBECCA_mod <- '
data {
int<lower=1> n;
int<lower=1> p;
int<lower=1> q;
int<lower=1> k;
int<lower=1> k1;
int<lower=1> k2;
int<lower=1> D;
int<lower=1> K;
int<lower=0, upper=2> X[n,p];
real Y[n,q];
vector[K] mu_z;
cov_matrix[K] Sigma_z;
#real<lower=0> lambda;
}
parameters {
real<lower=0> alpha[K]; # top part of E
real<lower=0> beta[K]; # bottom part of E
real<lower=0> sigma[q]; # modeling data Y, Y is normally distributed
matrix[n,K] Z;
real<lower=0> matrix[D,k] Pi;
matrix[D,K] H;
matrix[D,K] E;

}
transformed parameters {
matrix[D,K] Phii;
matrix[n,D] Theta;
matrix[n,p] Theta_x;
matrix[n,q] Theta_y;
matrix[n,p] p_x;
matrix[n,q] mu_y;

Phii <- E .* H;
Theta <- Z*transpose(Phii);

for(i in 1:n)
{for(j in 1:p)
Theta_x[i,j] <- Theta[i,j];

for(l in 1:q)
Theta_y[i,l] <- Theta[i,l+p];}

{for(i in 1:n)
for(j in 1:p)
p_x[i,j] <- exp(Theta_x[i,j])/(1+exp(Theta_x[i,j]));}

{for(i in 1:n)
for(l in 1:q)
mu_y[i,l] <- Theta_y[i,l];}

}
model {
############## Prior #################

for Z

for(i in 1:n)
{Z[i] ~ multi_normal(mu_z, Sigma_z);}

for alpha, beta (E)

for(j in 1:K)
{alpha[j] ~ inv_gamma(0.01, 0.01);
beta[j] ~ inv_gamma(0.1, 0.1);} ## This is just an initial value

for H

for(i in 1:D)
for(j in 1:K)
Pi[i,j] ~ exponential(2);
for(i in 1:D)
for(j in 1:K)
H[i,j] ~ normal(0,Pi[i,j]);

for(j in 1:K)
for(i in 1:p)
E[i,j] ~ normal(0, alpha[j]);

for(j in 1:K)
for(l in 1:q)
E[l,j] ~ normal(0, beta[j]);

for variance of y

for(j in 1:q)
sigma[j] ~ inv_gamma(0.01, 0.01);

Generate Data

for(i in 1:n)
{
for(j in 1:p)
{
X[i,j] ~ binomial(2,p_x[i,j]);
}
for(l in 1:q)
{
Y[i,l] ~ normal(mu_y[i,l], sigma[l]);
}
}
}
'
SBECCA_dat<-list(“X”=X,“Y”=Y, “mu_z” = rep(0,7), “Sigma_z” = diag(rep(1,7)),“n”=100,“p”=10,“q”=9,“k”=4, “k1”=2, “k2”=1, “D”= 19, “K”=7)
happyfit <- stan(model_code = SBECCA_mod, data = SBECCA_dat, thin=3, iter = 100, chains =2)

The variable Pi is defined as:

real<lower=0> matrix[D,k] Pi;

But the loop is over:

for(i in 1:D)
  for(j in 1:K)
    Pi[i,j] ~ exponential(2);
  }
}

Is the difference between k and K the problem (Stan is case sensitive)? Also, if you surround your code on Discourse with a pair of triple backticks and you can get it to format nicely. It’s hard to read without that.

Since you’re talking about horseshoe, is it that you wanted to use some sort of Bayesian Lasso/Laplace/double_exponential prior (if this is for Pi, there shouldn’t be a lower bound of zero on it)? @anon75146577 made a blog post about that: The king must die | Statistical Modeling, Causal Inference, and Social Science . Maybe you could find some answers there.