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)