Create augmented_ilr.stan

[spinkney]

As discussed in #39 and #43

[spinkney]

Should I also code up the Vinv matrix in Stan? Or can we leave that to be calculated outside? @mjhajharia @bob-carpenter

[bob-carpenter]

If it's something needed for the transform, I'd prefer to see it calculated in the Stan program. We don't have ways to pass data to transforms in the current version of Stan. Of course, if it's just a constant calculated once, we can and should do that on the outside and hard code it.

[spinkney]

@bob-carpenter @mjhajharia I think this can be included in the paper. It's simply the alr but with the Vinv which only affects the jacobian by the det of Vinv which is constant. Anyway, it would be interesting to see if it's "better" than alr.

[bob-carpenter]

If it works out to be the same, we can drop it out of figures and just include a footnote. Adding a constant won't change the curvature or derivatives---is there a chance it makes arithmetic more stable or does it center inits or something better?

[spinkney]

The ESS is doubled on the $N^{th}$ value, where I'm using the softmax.stan and the code in this pr with

sample(
  data = list(N = N),
  parallel_chains = 4,
  seed = 123124,
  iter_sampling = 500,
  iter_warmup = 500
)

> out_ilr
# A tibble: 5 × 10
  variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
  <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
1 x[1]     0.201  0.163 0.160 0.155 0.0143 0.518 1.00     1421.     806.
2 x[2]     0.203  0.169 0.159 0.160 0.0113 0.516 1.00     1742.    1135.
3 x[3]     0.205  0.166 0.168 0.161 0.0126 0.538 1.00     1690.    1257.
4 x[4]     0.201  0.158 0.168 0.156 0.0105 0.542 1.00     1779.    1059.
5 x[5]     0.190  0.150 0.157 0.146 0.0113 0.507 0.999    1249.     685.
> out_alr <- mod_alr_out$summary("x")
> out_alr
# A tibble: 5 × 10
  variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
  <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
1 x[1]     0.198  0.160 0.159 0.149 0.0128 0.528  1.00    1854.    1177.
2 x[2]     0.201  0.164 0.160 0.154 0.0111 0.527  1.00    1909.     998.
3 x[3]     0.199  0.160 0.163 0.161 0.0109 0.528  1.00    1639.    1011.
4 x[4]     0.197  0.159 0.161 0.152 0.0125 0.530  1.00    1611.    1008.
5 x[5]     0.205  0.171 0.163 0.159 0.0135 0.525  1.00     608.     581.

with 2500 samples

> out_ilr
# A tibble: 5 × 10
  variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
  <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
1 x[1]     0.198  0.156 0.163 0.151 0.0125 0.528  1.00    7929.    4761.
2 x[2]     0.199  0.159 0.163 0.156 0.0121 0.529  1.00    7354.    4072.
3 x[3]     0.202  0.164 0.163 0.157 0.0133 0.530  1.00    8472.    5015.
4 x[4]     0.200  0.162 0.163 0.154 0.0131 0.522  1.00    8435.    5417.
5 x[5]     0.200  0.160 0.163 0.155 0.0126 0.531  1.00    8669.    4825.
> out_alr
# A tibble: 5 × 10
  variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
  <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
1 x[1]     0.201  0.159 0.165 0.155 0.0122 0.530  1.00    7639.    4817.
2 x[2]     0.200  0.159 0.165 0.155 0.0125 0.535  1.00    8278.    5214.
3 x[3]     0.198  0.159 0.163 0.154 0.0121 0.524  1.00    7309.    5206.
4 x[4]     0.200  0.158 0.163 0.150 0.0132 0.525  1.00    7949.    4974.
5 x[5]     0.201  0.160 0.163 0.150 0.0132 0.527  1.00    3819.    2876.

[mjhajharia]

I haven't checked the code yet, but we use additive log ratio transform or ALR in the paper? was ILR a typo? they're different transforms, ILR is the isomorphic one

[mjhajharia]

is it the one mentioned in sec 5.4.2 here https://people.stat.sc.edu/hansont/stat730/paketo-libre.pdf