# Efficiently sampling a matrix of independent parameters

**URL:** <https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175>\
**Category:** Modeling\
**Tags:** performance\
**Created:** [May 12, 2018, 2:11am UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175 "2018-05-12T02:11:15Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![jayrobwilliams](https://avatars.discourse-cdn.com/v4/letter/j/91b2a8/32.png) [@jayrobwilliams](https://discourse.mc-stan.org/u/jayrobwilliams)\
**Post date:** [May 12, 2018, 2:11am UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175/1 "2018-05-12T02:11:15Z")

</div>

What is the most efficient way to sample a matrix of parameters drawn independently from the same distribution? I’m working on a multinomial logit model where I have a matrix of coefficients \beta that includes one column for each of K-1 outcomes in the data. I have \beta declared as a D \times K-1 matrix where each element is distributed \mathcal{N}(\mu\_\beta, \sigma\_\beta), where \mu\_\beta and \sigma\_\beta are hyperpriors.

I’m currently using `to_vector` to draw \beta because you can’t have a matrix on the left hand side of a `normal_lpdf`. A MWE of the model is included below. Is this the most efficient approach to sampling these parameters?

```stan
data {
  
  int<lower=1> N; // number of observations
  int<lower=1> D; // number of predictors (p + 1)
  int<lower=2> K; // number of alternatives
  int<lower=1,upper=K> Y[N]; // response variable
  vector[D] X[N]; // design matrix

}

parameters {
  
  matrix[K-1,D] beta_free_raw; // regression coefficeients for K - 1 alternatives
  real mu_beta;
  real<lower=.001> sigma_beta;
  
}

transformed parameters {
  
  matrix[K-1,D] beta_free;
  matrix[K,D] beta;
  beta_free = mu_beta + sigma_beta * beta_free_raw; // re-center parameters
  beta = append_row(rep_row_vector(0, D), beta_free); // append row of 0s to top
  
}

model {
  
  mu_beta ~ normal(0, 5);
  sigma_beta ~ cauchy(0, 2.5);
  to_vector(beta_free_raw) ~ normal(0, 1);
  
  for (n in 1:N) Y[n] ~ categorical_logit(to_vector(beta * X[n]));
  
}

```

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 12, 2018, 3:25am UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175/2 "2018-05-12T03:25:46Z")

</div>

Using `to_vector(beta_free_raw) ~ normal(0, 1);` but you should not read that as `beta_free_raw` is drawn from a standard normal distribution. Rather, the standard normal PDF in log form is being tacked on to the kernel for each element of `beta_free_raw`. That is why it is clearer to write `target += normal_lpdf(to_vector(beta_free_raw) | 0, 1);`.

---

<div class="post-metadata">

**Author:** ![jayrobwilliams](https://avatars.discourse-cdn.com/v4/letter/j/91b2a8/32.png) [@jayrobwilliams](https://discourse.mc-stan.org/u/jayrobwilliams)\
**Post date:** [May 13, 2018, 3:16pm UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175/3 "2018-05-13T15:16:25Z")

</div>

Got it. What are the practical differences between those two statements? Also, are there any other ways to write this model more efficiently?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 13, 2018, 4:09pm UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175/4 "2018-05-13T16:09:41Z")

</div>

The practical difference is that it does not confuse yourself or others as to what Stan is doing. That and the `target +=` version calculates constants, which is somewhat slower but necessary in some situations.

---

<div class="post-metadata">

**Author:** ![jayrobwilliams](https://avatars.discourse-cdn.com/v4/letter/j/91b2a8/32.png) [@jayrobwilliams](https://discourse.mc-stan.org/u/jayrobwilliams)\
**Post date:** [May 13, 2018, 4:42pm UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175/5 "2018-05-13T16:42:33Z")

</div>

Sorry if I wasn’t clear. I meant why is it important to remind yourself/others that you’re incrementing the log probability instead of drawing from a standard normal?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 13, 2018, 4:49pm UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175/6 "2018-05-13T16:49:17Z")

</div>

Because in order to sample from a posterior distribution in Stan, the Stan program has to commit to a log-kernel of a density function. And the shape of that log-kernel affects how well NUTS can draw from that posterior distribution. All of the tricks to improve the effective sample size are based on manipulating the log-kernel to get better geometry. When you ask “Also, are there any other ways to write this model more efficiently?” that is equivalent to asking “Also, are there any other ways to manipulate the log-kernel to get better geometry?”

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [June 15, 2018, 10:39pm UTC](https://discourse.mc-stan.org/t/efficiently-sampling-a-matrix-of-independent-parameters/4175/7 "2018-06-15T22:39:26Z")

</div>

> [@jayrobwilliams](#):
>
> why is it important to remind yourself/others that you’re incrementing the log probability instead of drawing from a standard normal?

It’s often possible (and useful) to think of a model generatively, but that’s not what the computation is doing. The computation is just computing a log density and its gradient and the geometry of that log density is what will determine statistical sampling efficiency.

> [@bgoodri](#):
>
> When you ask “Also, are there any other ways to write this model more efficiently?” that is equivalent to asking “Also, are there any other ways to manipulate the log-kernel to get better geometry?”

I’d qualify that. There are ways to make the iterations go faster (usually precomputation and vectorization) and there are ways to make the geometry better so that mixing is better (usually with reparameterizations).