E_loo function in loo package, type = "var"

General
1

I’m using importance sampling, via the loo package, to correct for an approximation of the likelihood. Using log ratios, I get (hopefully) improved Monte Carlo estimators, using

log_ratios <- fit$draws("log_ratios")
psis_fit <- psis(log_ratios, r_eff =  relative_eff(log_ratios))
mean <- E_loo(parm_draw, psis_fit, log_ratios = log_ratios, type = "mean")$value

In addition to the mean, I compute sample estimators of the variance and the quantiles. In my particular example, I trust the approximate likelihood and expect the original Monte Carlo estimators and the IS estimators to be in good agreement. This is true in all cases, except for the variance.

From Stan’s MCMC:

> fit$summary()
# A tibble: 582 × 10
   variable    mean  median      sd     mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <num>   <num>   <num>   <num>   <num>   <num> <num>    <num>    <num>
 1 lp__       9.04    9.41   1.89    1.71     5.43   11.4    1.00    1427.    2019.
 2 CL        16.6    16.6    0.444   0.443   15.9    17.3    1.00    3123.    2512.
 3 Q         23.0    22.9    1.82    1.77    20.1    26.1    1.00    2788.    2270.
 4 VC        34.1    34.5    5.10    4.65    24.6    41.8    1.01    1577.    1231.
 5 VP       235.    234.    20.7    20.0    204.    272.     1.00    3101.    2287.
 6 ka         2.77    2.77   0.547   0.513    1.83    3.67   1.01    1576.    1152.

From the IS estimators:

  parms        mean          var         q5       q95         khat
1    CL  16.6282713 5.621361e-05  15.895600  17.33510  0.038601727
2     Q  22.9959122 9.481787e-04  20.074000  26.14422 -0.132261381
3    VC  34.0892628 7.432019e-03  24.605302  41.80701 -0.068656790
4    VP 235.4317697 1.219584e-01 203.626003 271.55000 -0.040339663
5    ka   2.7660148 8.533706e-05   1.823910   3.67362  0.102390467
6 sigma   0.1675363 9.074186e-08   0.140867   0.19896  0.007235841

Based on the orders of magnitude I get for var this seems to be an estimate of the MCSE (although it doesn’t quite match the MCSE returned by Stan…).

Am I misunderstanding type = "var"?

2 Likes
2

Can you be more specific? cmdstanr? rstan?

Whatever is going on, the column title should be changed. This is obviously not the posterior variance of the parameters given the 90% posterior intervals.

Also, we should always be reporting on the sd scale, not the var scale.

The good news is that the IS didn’t really change anything as far as I can see, which makes sense given the very low k-hat estimates.

3

Hmmm, let’s check with @avehtari.

Also, here’s are internal functions used to calculate this (w are the IS weights, not on log scale). Nothing jumps out at me as obviously wrong, but I’m a bit pressed for time at the moment so I could be overlooking a bug here (or in the code that wraps these functions):

4

I couldn’t quite follow that code, but given it’s dividing by r_eff, it’s not variance. From the code, it looks like an importance sampled estimate of squared standard error.

P.S. Calculating variance with x - mean(x) can be unstable due to catastrophic cancellation if the value of the variable is much larger than its standard deviation. Having said that, I think Stan’s calculations suffer from the same problem and it’s only a problem with edge cases.

5

I agree that’s what it looks like. It’s been years now so I actually forget the original intent, but hopefully @avehtari will remember and can tell us if it’s just a naming issue or if we’re not computing what we originally intended.

I hadn’t thought about that before but yeah that makes sense!

6

Oops. It’s no computing what we originally intended. It’s not just the r_eff part (which by default has no effect as by default it’s 1), but the sum(w^2 * r) part is clearly from the variance for the mean part. .wvar should return

(sum(.wmean(x^2, w))-sum(.wmean(x, w)^2))/(1-1/length(w))

Which would match base R var() in case of equal weights.

We could also add support for type="sd" which would be just a square root of that.

@charlesm93 thanks for reporting this

2 Likes
7

Just opened a PR to fix this:

2 Likes
8

@charlesm93 I think you said you wanted to use this in a workshop and I’m not sure when we’ll submit this to CRAN (hopefully soon), but the fix is now merged and available by installing the master branch from GitHub: remotes::install_github("stan-dev/loo"). So if we don’t get it on CRAN in time for your class you can have them install from GitHub with that one line of code.

1 Like
9

And we added type = "sd" as an option.

1 Like