I am fitting a model where unobserved parameters in a 2-simplex are transformed from a multivariate normal, then they generate the actual data. The relevant lines are not unlike

```
for (i in 1:N) {
z[i] ~ multi_normal_cholesky(mu, L); // 2-dimensional mv normal
a[i] = inv_logit(z[i][1]); // transform to [0,1]
b[i] = (1-a[i])*inv_logit(z[i][2]); // transform to make a simplex
// ... a[i] and b[i] unobserved, they are parameters to a multinomial
// ... that generates observations, omitted here
}
```

The multivariate normal is just a modeling device, I care about `a`

and `b`

. `N`

is around 1000, so the questions is how to visualize this. I thought of using HPD regions graphically. I could

- use the posterior mean of
`mu`

and`L`

, and visualize the conditional`a`

and`b`

based on that, or - I could use
*all*the`a`

and`b`

s, and plot HPD regions for that.

I am wondering which one makes more sense conceptually. (1) does not capture posterior uncertainty, while I think that (2) would overstate it.