Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] draft subsampling bootstrap for mcse #1974

Draft
wants to merge 6 commits into
base: main
Choose a base branch
from

Conversation

OriolAbril
Copy link
Member

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

  • Follows official PR format
  • Includes a sample plot to visually illustrate the changes (only for plot-related functions)
  • New features are properly documented (with an example if appropriate)?
  • Includes new or updated tests to cover the new feature
  • Code style correct (follows pylint and black guidelines)
  • Changes are listed in changelog

@OriolAbril OriolAbril marked this pull request as draft February 3, 2022 13:32
https://doi.org/10.1214/14-EJS957

"""
flat_ary = np.ravel(ary)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this sensitive for which order the ravel is done?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@OriolAbril have you tested alternative approaches for handling multiple chains yet?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I benchmarked 3 different approaches for estimating the mcse of the mean

  • ess: using ess to estimate mcse, for reference
  • sbm_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 by nchains^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:
    mcse_methods_rmse

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.

arviz/stats/diagnostics.py Outdated Show resolved Hide resolved
@OriolAbril
Copy link
Member Author

This func method for mcse should also replace mc_error

Comment on lines 418 to 422
if prob is not None:
func_kwargs["prob"] = prob
elif func is not None:
func_kwargs["func"] = func

Copy link
Contributor

@ahartikainen ahartikainen Mar 14, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

setdefault ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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))
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link

codecov bot commented Aug 17, 2022

Codecov Report

Merging #1974 (13dd989) into main (8a2bc39) will decrease coverage by 0.16%.
The diff coverage is 36.84%.

@@            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     
Impacted Files Coverage Δ
arviz/stats/diagnostics.py 93.28% <36.84%> (-5.69%) ⬇️
arviz/data/datasets.py 98.48% <0.00%> (+0.09%) ⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants