I’d like to inquire about how to accelerate the sampling process in Stan for my multi-layer hierarchical model. Let me briefly describe it: In the first layer, I draw n values from a beta distribution. In the second layer, I perform a variable transformation to convert these n variables into (alpha, beta) pairs. Then, I apply the betabinomial model (alpha_n, beta_n, times) to these transformed values. In simple terms, I want to use sampling to find the hierarchical set of n parameters. However, currently, when simulating the model, I can only use 2.5% of n (500/20000) to make it computationally feasible. I’d appreciate any advice on methods to speed up the process. Below is the code I’m using. Thank you all for your help.

The math form of my model:

\mu_i\sim Beta(\alpha_u, \beta_u)

trial_i\sim Gamma(10,10),

where pdf of gamma is

f(x|\alpha, \beta)=\frac{1}{Γ(\alpha)\beta^\alpha}x^{1-\alpha}e^{-\frac{x}{\beta}}

With parameter tranformation:

\alpha_i=\mu_i*trial_i

\beta_i=(1-\mu_i)*trial_i

The relationship between datas and parameters:

ctr_{ij}\sim Beta(\alpha_i, \beta_i)

click_{ij}\sim Bin(impression_{ij}, ctr_{ij})

Let i denote the i’th merchant, let (i,j) denote the j’th kind of goods from i’th merchant

The code is here

```
data{
int<lower=0> N;
int<lower=0> M;
int<lower=1, upper=466> type[N];
int<lower=0> impression[N];
int<lower=0> click[N];
}
parameters{
vector<lower=0.01, upper=1>[M] mu;
vector<lower=0>[M] trial;
real<lower=0> alpha_u;
real<lower=0> beta_u;
}
transformed parameters{
vector[M] alpha_i;
vector[M] beta_i;
alpha_i = mu.*trial;
beta_i = (1-mu).*trial;
}
model{
alpha_u~uniform(0.01, 200);
beta_u~uniform(0.01, 200);
mu~beta(alpha_u, beta_u);
trial~gamma(10, 0.1);
for(i in 1:N) {
click[i]~beta_binomial(impression[i], alpha_i[type[i]], beta_i[type[i]]);
}
}
```

Besides, these are the parameters.

```
stan_data = {'N':int(ruten_sample.shape[0]),
'M':int(ruten_sample.seller_number.nunique()),
'type': np.array(ruten_sample['seller_number']).astype(int),
'impression': np.array(ruten_sample['impression']).astype(int),
'click': np.array(ruten_sample['click']).astype(int)}
posterior = stan.build(text, data=stan_data)
fit = posterior.sample(num_chains=1, num_samples=2000)
df = fit.to_frame()
```