I would like to implement sparse infinite factor priors (ala Bhattacharya & Dunson 2011) on the loading matrix in a latent factor model. Sparse infinite factor priors allow for introducing infinitely many factors, with the loadings increasingly shrunk towards zero as the column index increases. I’m very new to Stan, but implemented this in Jags as follows

```
For j=1,...,J rows and q=1,...,Q latent factors for loading matrix L
for ( j in 1:J ) { for ( q in 1:Q ) { L[j,q] ~ dnorm( 0, phi[j,q]*tau[q] ) } }
for ( j in 1:J ) { for ( q in 1:Q ) { phi[j,q] ~ dgamma( 3/2, 3/2 ) } }
for ( q in 1:Q ) { tau[q] = prod( delta[1:q] ) }
for ( l in 1:Q ) { delta[l] ~ dgamma( 50, 1 ) }
(rate and shape parameters taken from Ovaskainen et al. 2015)
```

My poor attempt to do the same in Stan

```
// Toy model
data {
int<lower=1> N; // number of observations
int<lower=1> J; // dimension of observed data (i.e., number of cols in y_ij)
int<lower=1> I; // number of sample units (i.e., number of rows in y_ij)
int counts[N]; // counts in y_ij
int<lower=1> Q; // number of latent factors
int <lower=1,upper=I> row_id[N]; // row indices for observations
int <lower=1,upper=J> col_id[N]; // column indices for observations
}
parameters {
matrix[I,Q] LF;
matrix[J,Q] L;
matrix[J,Q] phi;
vector[Q] delta;
}
transformed parameters {
real log_lambda[N];
vector[Q] tau;
for ( q in 1:Q ) tau[q] = prod(delta[1:q]);
//linpred
for ( n in 1:N ) log_lambda[n] = LF[row_id[n],]*L[col_id[n],]';
}
model {
for( i in 1:I ) for ( q in 1:Q ) LF[i,q] ~ normal(0,1);
for ( j in 1:J ) for ( q in 1:Q ) L[j,q] ~ normal( 0, phi[j,q]*tau[q] );
for ( j in 1:J ) for ( q in 1:Q ) phi[j,q] ~ gamma( 3/2, 3/2 );
for ( l in 1:Q ) delta[l] ~ gamma( 50, 1 );
counts ~ poisson_log( log_lambda );
}
```

This produces 100 error msgs:

```
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Scale parameter is -0.00111817, but must be > 0!
Initialization between (-0.1, 0.1) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
```

Data-generation model in R

```
J=20
I=100
Q=5
log_lambda <- y <- matrix(0,nrow=I,ncol=J)
beta0 <- runif(J,-1,1)
alpha0 <- runif(I,-1,1)
#simulating sparse factor priors
LF <- matrix(0,nrow=I,ncol=Q)
for(k in 1:2) {LF[,k] <- rmvnorm(1, mean = rep(0,I))}
L <- matrix(0,nrow=J,ncol=Q)
for(k in 1:2) {L[,k] <- rmvnorm(1, mean = rep(0,J))}
for(i in 1:I) {
for(j in 1:J) {
for(q in 1:Q) {
log_lambda[i,j] <- exp(LF[i,q]%*%t(L[j,q))
y[i,j] <- rpois(1,log_lambda[i,j])
}}}
#prepare data for stan
yy <- reshape2::melt(y)
colnames(yy) <- c("rows","cols","counts")
row_id <- yy$rows
col_id <- yy$cols
counts <- yy$counts
#run stan
stan_d <- list(I = I, J = J, N = N, Q = Q, counts = counts, row_id = row_id, col_id = col_id)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
m_mod <- stan_model("model.stan")
m_fit <- sampling(object = m_mod, data = stan_d, chains = 1, iter = 2000,
pars = c(""log_lambda","L","LV","phi","tau","delta"), init_r = .1)
```

Any advice is highly appreciated.