[WIP] draft subsampling bootstrap for mcse #1974
Conversation
| https://doi.org/10.1214/14-EJS957 | ||
|
|
||
| """ | ||
| flat_ary = np.ravel(ary) |
There was a problem hiding this comment.
Is this sensitive for which order the ravel is done?
There was a problem hiding this comment.
I think it technically is, but there should be no difference (hopefully) if the model has converged. It is also not clear to me how should multiple chains be handled when implementing this algorithm, I started with this flatten approach but I can test a couple options.
There was a problem hiding this comment.
@OriolAbril have you tested alternative approaches for handling multiple chains yet?
There was a problem hiding this comment.
I benchmarked 3 different approaches for estimating the mcse of the mean
ess: usingessto estimate mcse, for referencesbm_stack_chains: estimating mcse with SBM by concatenating the chains (as done here)sbm_stack_draws: estimating mcse with SBM by interleaving chains (i.e. flattening on the other dimension)sbm_separate_chains: estimate the mcse of each chain with SBM and then sum the variances and divide bynchains^2.sbm_shuffle: flatten the chains and shuffle the draws before running SBM
I performed the same benchmark as in https://avehtari.github.io/rhat_ess/ess_comparison.html, and transforming the chains to target different stationary distributions. Here's the result:
In all cases SBM underestimates the MCSE; this is particularly severe when autocorrelation is high and sample sizes are low. sbm_stack_chains is consistently better than the alternatives though. I didn't even bother plotting sbm_shuffle, since it was apparent pretty quickly that it was far worse than the others.
|
This |
arviz/stats/diagnostics.py
Outdated
| if prob is not None: | ||
| func_kwargs["prob"] = prob | ||
| elif func is not None: | ||
| func_kwargs["func"] = func | ||
|
|
There was a problem hiding this comment.
nevermind, bad suggestion
| for i in range(n - b): | ||
| sub_ary = flat_ary[i : i + b] | ||
| func_estimates[i] = func(sub_ary, **func_kwargs) | ||
| func_estimate_sd = np.sqrt(b * var_func(func_estimates)) |
There was a problem hiding this comment.
we should probably decide API-wise if we want to keep this or instead move to std_func and multiply that by the square root of b.
| Quantile information. | ||
| func : callable, optional |
There was a problem hiding this comment.
we could also consider allowing some strings here. e.g. using "circmean" expands to stats.circmean as func here and also fills the mcse_kwargs with {"var_func": stats.circvar}
Codecov Report
@@ Coverage Diff @@
## main #1974 +/- ##
==========================================
- Coverage 90.78% 90.62% -0.17%
==========================================
Files 117 117
Lines 12484 12518 +34
==========================================
+ Hits 11334 11344 +10
- Misses 1150 1174 +24
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. |
Description
cc @aloctavodia
This adds a draft version of the subsampling bootstrap method (as defined in https://doi.org/10.1214/14-EJS957) but allowing for it to be used on any arbitrary function.
Checklist