I am currently learning stan from a book called “Bayesian Models for Astrophysical Data” and I am totally hooked. However when I ran stan on a simple linear mixed model

y_{i,j} = \beta_0 + \mathbf{x}_{i,j}^T\beta + \alpha_i + \epsilon_{i,j}

with 4500 observations and a total of 3 fixed effect and 20 random effect coefficients (+ 2 coefficients from a prior on each group of coefficients) it took quite some time (~4.7mins). I used the following mcmc config

`fit = pystan.stan(model_code=stan_code, data=model_data, iter=10000, chains=3, thin=10, warmup=6000, n_jobs=3)`

In one of the applications where I want to use the model I have 1.1m observations and probably at least 10k parameters (because there are 2 hierarchies with 3000 and 1000 categories each). Is this doable at all in stan?

Thanks!