Hi stan community,

I am trying to run a hierarchical survival model in Stan, but I found sampling is very time-consuming (takes more than 10 days even for 25000 iteration/10000 warmup).

I am not good at computing, so would be greatly appreciated if there are any tips!

My stan code is as below (N=34810, n_grid=350, province=31, num=12).

```
data {
int <lower=0> province;
int <lower=0> n_grid;
int <lower=0> N;
int <lower=0> num;
int <lower=0> number;
vector <lower=0> [province] p;
matrix [province,n_grid] kernel;
matrix [N,num] predictor;
matrix [N,number] respon;
}
parameters{
vector [num] beta;
real <lower=0> sigma0;
real <lower=0> sigma;
vector [n_grid] x;
}
transformed parameters{
real alpha = 1 / sigma;
vector [N] lambdaa;
vector[province] a;
for(n in 1:N) {
lambdaa [n] = exp ((-1 * predictor [n] * beta) / sigma);
}
{ // `z` is only visible inside these braces, not part of output
vector[province] z;
for (k in 1 : province) {
vector[n_grid] landa_k;
for (j in 1 : n_grid) {
landa_k[j] = kernel[k, j] * exp(x[j]);
}
z[k] = sum(landa_k) * p[k];
}
a = z / sum(z); // assign `a`
}
}
model{
matrix [N,province] f;
matrix [N,province] s;
for (j in 1:n_grid){
target += normal_lpdf(x[j]|0,1);
}
for(u in 1:num){
target += normal_lpdf(beta [u]| 0, sigma0);
}
target += inv_gamma_lpdf(sigma0| 1, 1);
target += normal_lpdf(log (alpha) | 0, 5);
for (i in 1:N) {
for (k in 1:province){
f[i,k]= a[k]*((lambdaa [i]*alpha*respon[i,1]^(alpha-1))/((1+(lambdaa [i]*respon[i,1]^alpha))^2));
s[i,k]= a[k]*(1- ((lambdaa [i]*respon[i,1]^alpha)/ (1+(lambdaa [i]*respon[i,1]^alpha))));
target += (respon[i,2] * log (f[i,k])) + ((1 - respon[i,2])* log(s[i, k]));
}
}
}
```

Any tips and comments are welcome. Thanks in advance!

Best,

Eisa

**moderator** edited for stan code formatting