# Optimize Stan code - hierarchical model

**URL:** <https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646>\
**Category:** Modeling\
**Created:** [March 11, 2020, 9:31am UTC](https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646 "2020-03-11T09:31:28Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [March 11, 2020, 9:31am UTC](https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646/1 "2020-03-11T09:31:28Z")

</div>

Hi,

Is there any better way to parametrize this hierarchical model?

Thank you,

```stan
data {
 int<lower=0> N; //number of observations
 real<lower=0> CPS[N]; //Cost per service (response variable) 
 int<lower=0> Njkm; //total number of state-service-market (168 total, 28 do not have observed CPS)
 int<lower=0> Nkm; //number of market services (8)
 //int<lower=1> jkm_identifier[N]; //state-service-market identifier level 1
 //int<lower=1> km_identifier[N]; //service-market identifier level 1
 //int<lower=0> gender_identifier[N]; //population gender effect
 matrix[N,Njkm] design_matrix_1;
 matrix[N,Nkm] design_matrix_2;
 matrix[Njkm,Nkm] design_matrix_3;
 matrix[Njkm,24] design_matrix_4;
}
transformed data {
  }
parameters {
// real beta_0;
 vector[Njkm] mu_jkm;
 vector[Nkm] gamma_km;
 vector[24] beta_km;
 vector<lower=0>[Nkm] sigma_cps_km;
 vector<lower=0>[Nkm] sigma_mu_km;
 //real theta_g; //Population gender effect
}
transformed parameters {
     vector[Njkm] gamma_km_index = design_matrix_3 * gamma_km + design_matrix_4 * beta_km ; 
     vector<lower=0>[Njkm] sigma_mu_index = design_matrix_3 * sigma_mu_km;
     
     vector[N] mu_index =design_matrix_1 * mu_jkm;
     vector<lower=0>[N] sigma_cps_index = design_matrix_2 * sigma_cps_km;

    for (n in 1:N){
      mu_index[n]=exp(mu_index[n]);
    }
}
model { 
//Priors
//target += normal_lpdf(beta_0jkm | beta_km, sigma_u0jkm);

//likelihood
//for (i in 1:N){
//CPS[i] ~ normal (mu[i], sigma_cps[i]);
target += normal_lpdf(CPS |mu_index, sigma_cps_index);
target += normal_lpdf(mu_jkm |gamma_km_index, sigma_mu_index);
}
generated quantities{
   real yrep[N];
   for (i in 1:N)
   yrep[i]= normal_rng(mu_index[i], sigma_cps_index[i]);
}
```

---

<div class="post-metadata">

**Author:** ![emiruz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/emiruz/32/5518_2.png) [@emiruz](https://discourse.mc-stan.org/u/emiruz)\
**Post date:** [March 11, 2020, 5:28pm UTC](https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646/2 "2020-03-11T17:28:12Z")

</div>

So in [this thread](https://discourse.mc-stan.org/t/increasing-performance-of-this-code-for-logistic-regression/13389)@avehtari gives OP the following advise:

> [@Increasing performance of this code for logistic regression](https://discourse.mc-stan.org/t/increasing-performance-of-this-code-for-logistic-regression/13389/3):
>
> You have not set any prior on coefs, which means the prior is improper and combined with bernoulli\_logit the posterior can be improper and MCMC can’t produce valid posterior draws even in infinite time. Add proper priors and report the result here

if you’re looking to reparameterise because the model runs slow or doesn’t converge then @avehtari advise RE priors might apply here too since none of the parameters seem to have priors in your model unless I’ve misunderstood something.

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [March 11, 2020, 5:34pm UTC](https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646/3 "2020-03-11T17:34:20Z")

</div>

Yes, thank you.  
Another point, is this matrix notation the most efficient way to write the problem?

---

<div class="post-metadata">

**Author:** ![emiruz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/emiruz/32/5518_2.png) [@emiruz](https://discourse.mc-stan.org/u/emiruz)\
**Post date:** [March 11, 2020, 5:45pm UTC](https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646/4 "2020-03-11T17:45:31Z")

</div>

Just IMHO, and assuming your objective is performance, I would tend to think about it in terms of (1) Am I performing any operations that I don’t need to? (2) Have I vectorised as much as possible?

According to [this](https://mc-stan.org/docs/2_20/functions-reference/fun-vectorization.html), the exp function is vectorised. so you can write :

```stan
mu_index=exp(mu_index);

```

instead of:

```stan
for (n in 1:N){
      mu_index[n]=exp(mu_index[n]);
    }

```

Similarly, [this](https://mc-stan.org/docs/2_20/functions-reference/normal-distribution.html) says that normal\_rng is vectorised so you can write this:

```stan
yrep=normal_rng(mu_index, sigma_cps_index);

```

instead of this:

```stan
 for (i in 1:N)
   yrep[i]= normal_rng(mu_index[i], sigma_cps_index[i]);

```

Except for that, I don’t see anything else personally, but others might.

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [March 11, 2020, 9:25pm UTC](https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646/5 "2020-03-11T21:25:19Z")

</div>

I can’t vectorize yrep:

## SYNTAX ERROR, MESSAGE(S) FROM PARSER: Variable definition base type mismatch, variable declared as base type vector variable definition has base type real[] error in ‘model49c4729690b\_stan\_CPS\_11’ at line 50, column 54

```
48: }
49: generated quantities{
50: vector[N] yrep = normal_rng(mu_index, sigma_cps_index);
                                                         ^
51: }

```

* * *

Error in stanc(file = file, model\_code = model\_code, model\_name = model\_name, :  
failed to parse Stan model ‘stan\_CPS\_11’ due to the above error.

---

<div class="post-metadata">

**Author:** ![emiruz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/emiruz/32/5518_2.png) [@emiruz](https://discourse.mc-stan.org/u/emiruz)\
**Post date:** [March 11, 2020, 9:31pm UTC](https://discourse.mc-stan.org/t/optimize-stan-code-hierarchical-model/13646/6 "2020-03-11T21:31:12Z")

</div>

Just make it an array instead ?