Computing Monte Carlo Standard Error for an arbitrary function

[MansMeg]

Hi!

Me, @paul-buerkner and @avehtari discussed the need for computing the MCSE for different functions of interest, such as the marginal variances, skewness and kurtosis of individual parameters. I think this functionality is probably best to implement directly into the posterior package. Aki suggested to use the batch means and subsampling approaches of https://projecteuclid.org/download/pdfview_1/euclid.ejs/1417615759 and I could implement this into posterior since I need to implement it anyway.

The most important part for me to know is how you would like the API to look. This is my first suggestion, feel free to comment or come up with a better suggestion.

monte_carlo_standard_error(x, FUN = mean)

where x is a draws object and FUN is an arbitrary function that operates on the individual parameters (rv?) of the draws object. It returns an object similar to summary() but with column with the names: [parameter_name]_[function_name] and [parameter_name]_[function_name]_se

Any thoughts or comments before I get started?

[paul-buerkner]

The interface should be similar to the existing mcse functions, which is that is only takes the draws of a single variable and returns a single value, the mcse of that variable for the given quantity. For the new function, we could use for instance mcse_fun(x, fun) as the function signature and also return a single value. This function could then be passed to summarise_draws to loop it over all quantities. For more inspriation see, for instance, ?mcse_mean

[MansMeg]

Great!

Then I guess

mcse_fun(x, FUN, method = "batch_means")

would be the way to go (I forgot the method arguments previously) and the function argument is using the R convention of FUN with capital letters (as in the *apply family)?

[paul-buerkner]

We rather follow tidyverse naming conventions than base R conventions so I would prefer to avoid any all capital argument names. We can debate on what would be a good argument name, but I would tend to oppose FUN.

[MansMeg]

Alright! Then I go with fun if thats the convention in posterior. I guess this is the styleguide? https://style.tidyverse.org/syntax.html

[paul-buerkner]

We didnt use an explicit style guide I think but more of our feeling which hopefully aligns with the style guide. :-D

I am open for other argument names of course just wanted to note that FUN may go against what we do elsewhere in posterior.

Måns Magnusson notifications@github.com schrieb am Fr., 3. Jan. 2020, 16:04:

Alright! Then I go with fun if thats the convention in posterior. I guess this is the styleguide? https://style.tidyverse.org/syntax.html

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jgabry/posterior/issues/57?email_source=notifications&email_token=ADCW2AGP7SZVZQQ4J24KQ7TQ35HZBA5CNFSM4KCNYTIKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEIBKAKI#issuecomment-570597417, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADCW2AHB4BHC3PZB6QFB5ALQ35HZBANCNFSM4KCNYTIA .

[mjskay]

I think getting this to work with summarise_draws as suggested will require allowing the measures argument to accept functions? (see #58). I like that approach.

Thinking about general use, this does create a bit of boilerplate in some ways, as you probably will have to give the output column a nice name:

summarise_draws(x, measures = c(mean, mean_se = ~ mcse_fun(., mean)))

I wonder if there is a cleaner way to support common patterns like that. E.g. a function like mcse_funs(funs) that takes a list of functions and returns a named list. So mcse_funs(c(mean, median)) returns:

list(
  mean_se = function(...) mcse_fun(fun = mean, ...),
  median_se = function(...) mcse_fun(fun = median, ...)
)

So you can do:

summarise_draws(x, measures = c(mean, median, mcse_funs(c(mean, median)))

Or something like that.

[MansMeg]

Thats great! Ill try to implement mcse_fun() asap so you can play around with it.

[MansMeg]

Alright. During the weekend I read up on MCSE and I think that the method we should use is the batch means (BM) methods with batch size \sqrt{S}, where S is the sample size. This is recommended below. The only problem is how to handle multiple chains in the best way.

I guess the simplest is to use the batch size b=sqrt(S_c) where S_c is the samples from chain c and then just combine all batches over all chains. I think this should not be problematic and the BM estimator should be consistent still for the MCMC SE for arbitrary (\pi integrable) function g.

Is this sufficient or should I try to dig in more?

See: Jones, Galin L., et al. "Fixed-width output analysis for Markov chain Monte Carlo." Journal of the American Statistical Association 101.476 (2006): 1537-1547.

[paul-buerkner]

Nice! Let's start as you suggested. Even if we change up the details with regard to the chains later on, most part of the code will remain the same anyway.

[MansMeg]

Yes, indeed.

Me and Aki spoke and he also wanted to check if we should make a smarter choice for batch sizes.

This is a small mini-project in itself, but I will refactor the Geyers http://github.com/jgabry/posterior/blob/master/R/convergence.R (lines 512-548) so it can be used for batch size choice.