Hi,

I am trying to fit the following hierarchical model where gene expression is estimated over many samples. Currently I have about 300 samples. When I run this with 12 genes for a test it runs about 7 minutes on our server that has a 32-core 2.4GHz CPU and 512GB RAM. But when I run it with full data (about 8,000 genes) it is taking a prohibitively long time (I had to kill it after running for a week).

I have already tried non-centered reparametrization for the lognormal part which was even a bit slower (with 12-gene test) than the one I am showing below. Could you help me scaling up this model? Any advice would be greatly appreciated.

```
data {
int<lower=1> M; // number of genes
int<lower=1> N; // number of samples
real<lower=0> C[N]; // size factor per sample
int<lower=0> X[M, N]; // observed count matrix
}
parameters {
real<lower=0> lambda[M, N];
real<lower=0> shape[M];
real<lower=0> scale[M];
real mu;
real<lower=0> sigma;
}
transformed parameters {
real<lower=0> rate[M];
real<lower=0> adj_lambda[M, N];
rate = inv(scale);
for (i in 1:M)
for (j in 1:N)
adj_lambda[i][j] = lambda[i][j] * C[j];
}
model {
for (i in 1:M) {
X[i] ~ poisson(adj_lambda[i]);
lambda[i] ~ gamma(shape[i], rate[i]);
}
shape ~ gamma(0.01, 0.01);
scale ~ lognormal(mu, sigma);
sigma ~ cauchy(0, 5);
}
```

I also wonder what you think about providing point estimators for `shape`

and `scale`

as in empirical Bayes perhaps with some large variance.

Thanks in advance.