Redesign of MCSE #63

sethaxen · 2023-01-16T12:53:46Z

Following #39 (comment) and #53, this PR redesigns mcse:

mcse now requires the user provide an estimator to make it less ambiguous what is estimated.
the vector methods to mcse have been removed and replaced with 3d array methods similar to ess_rhat, which also allow for eltypes to include Missing. This allows MCSEs to be computed for many parameters with one function call and also allows for improved ESS estimates when multiple chains are provided
a subsampling bootstrap method has been added as a fallback, since unlike batch means or overlapping batch means, it works for any estimator.
custom estimators have been added for mean, std, and quantile using the ESS methods
batch-means, initial monotone sequence, and initial positive sequence methods have been removed. These estimators only can estimate the MCSE of mean. From the Geyer reference in the docstrings, IMSE is better than IPSE, and both are better than batch-means, which tends to underestimate MCSE more severely. IMSE is used in the ess_rhat(::typeof(mean), x) call used in the new ess_rhat implementations.

Marked as a draft for now until I add tests.

coveralls · 2023-01-17T00:03:46Z

Pull Request Test Coverage Report for Build 3939266728

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

For more information on this, see Tracking coverage changes with pull request builds.
To avoid this issue with future PRs, see these Recommended CI Configurations.
For a quick fix, rebase this PR at GitHub. Your next report should be accurate.

Details

57 of 58 (98.28%) changed or added relevant lines in 5 files are covered.
No unchanged relevant lines lost coverage.
Overall coverage increased (+0.06%) to 95.828%

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
src/mcse.jl	39	40	97.5%

Totals
Change from base Build 3939418215:	0.06%
Covered Lines:	689
Relevant Lines:	719

💛 - Coveralls

codecov · 2023-01-17T00:04:54Z

Codecov Report

Base: 95.76% // Head: 96.16% // Increases project coverage by +0.40% 🎉

Coverage data is based on head (69fe0a1) compared to base (2e07d21).
Patch coverage: 100.00% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #63      +/-   ##
==========================================
+ Coverage   95.76%   96.16%   +0.40%     
==========================================
  Files          10       10              
  Lines         732      757      +25     
==========================================
+ Hits          701      728      +27     
+ Misses         31       29       -2

Impacted Files	Coverage Δ
src/ess.jl	`100.00% <100.00%> (ø)`
src/gewekediag.jl	`100.00% <100.00%> (ø)`
src/heideldiag.jl	`100.00% <100.00%> (ø)`
src/mcse.jl	`100.00% <100.00%> (+4.44%)`	⬆️
src/utils.jl	`100.00% <100.00%> (ø)`

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

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

sethaxen · 2023-01-17T00:33:40Z

src/mcse.jl

+    T = promote_type(eltype(x), typeof(zero(eltype(x)) / 1))
+    values = similar(x, T, (axes(x, 3),))
+    for (i, xi) in zip(eachindex(values), eachslice(x; dims=3))
+        values[i] = _mcse_sbm(f, vec(xi); batch_size=batch_size)


I haven't seen in the literature a study of how best to combine MCMC chains when estimating MCSE with SBM. I benchmarked a few alternatives and found this one underestimated the MCSE the least: arviz-devs/arviz#1974 (comment)

sethaxen · 2023-01-17T12:40:39Z

test/mcse.jl

-        for size in (51, 75, 100, 153)
-            @test_logs (:warn,) mcse(samples; method=:bm, size=size)
+@testset "mcse.jl" begin
+    @testset "estimand is within interval defined by MCSE estimate" begin


For mean, std, and quantile, these tests are somewhat redundant with those for ess_rhat, since the only difference is that for ess_rhat we compute the exact asymptotic variance of the estimator since we know the stationary distribution, whereas for mcse we estimate the asymptotic variance of the estimator. If we want to speed up our test suite, it may be better to remove the corresponding ess_rhat tests and use these as integration tests.

I've opted to speed up our test suite by decreasing the number of replicates (params) for each suite to 100 instead of 1,000.

sethaxen · 2023-01-17T12:44:53Z

src/mcse.jl

+    SBM tends to underestimate the MCSE, especially for highly autocorrelated chains.
+    SBM should only be used as a fallbeck when a specific [`mcse`](@ref) method for


This begs the question whether we should include mcse_sbm at all. On the one hand it's useful to have a fallback, and I don't know of any other general-purpose methods for estimating MCSE. On the other hand, we would prefer a fallback that errs on the side of overestimating MCSE instead of overestimating.

I lean towards including it with this caveat clearly documented.

Since we only want it to be used as a fallback, maybe we do not even want to add another function for it but also call it mcse? And just document there that it is used as a fallback but one should be aware of its limitations.

On second thought, it's useful to have it as a standalone method so if we want, we can compare its results with those of the specific estimators. But we can make the function internal and copy the documentation to mcse.

sethaxen · 2023-01-18T15:23:08Z

src/ess.jl

@@ -480,7 +482,7 @@ function _expectand_proxy(::typeof(StatsBase.mad), x)
    return _expectand_proxy(Statistics.median, x_folded)
 end
 function _expectand_proxy(f::Base.Fix2{typeof(Statistics.quantile),<:Real}, x)
-    y = similar(x, Bool)
+    y = similar(x)


Necessary because otherwise our transformation discards the type of x, so a Float32 eltype x will get a Float64 ESS/R-hat/MCSE.

improves the quality of the estimates and reduces random failures

src/heideldiag.jl

src/mcse.jl

devmotion · 2023-01-18T21:01:52Z

src/mcse.jl

+    SBM tends to underestimate the MCSE, especially for highly autocorrelated chains.
+    SBM should only be used as a fallbeck when a specific [`mcse`](@ref) method for


Since we only want it to be used as a fallback, maybe we do not even want to add another function for it but also call it mcse? And just document there that it is used as a fallback but one should be aware of its limitations.

src/mcse.jl

src/utils.jl

Co-authored-by: David Widmann <[email protected]>

sethaxen added 18 commits January 16, 2023 13:24

Add mcse_sbm

2d82222

Update description of estimator

e4b067b

Add specialized estimators for mean, std, and quantile

34e0221

Remove vector methods, defaulting to sbm

6839cd7

Update docstring

7dbda2d

Fix bugs

b93ce5b

Update docstrings

790a99b

Update docstring

9a1d3d5

Move helper functions to own file

c0d5a94

Rearrange tests

cf908af

Update mcse tests

93d121e

Export mcse_sbm

fe99356

Increment minor version number with DEV suffix

d12648b

Merge branch 'main' into mcseupdate

465afac

Increment docs and tests version numbers

7dac85e

Add additional citation

9407369

Merge branch 'main' into mcseupdate

f95a066

Update diagnostics to use new mcse

88b6c41

sethaxen added 2 commits January 17, 2023 01:05

Increase tolerance of mcse tests

899711e

Increase tolerance more

01b8dbc

sethaxen commented Jan 17, 2023

View reviewed changes

sethaxen marked this pull request as ready for review January 17, 2023 00:33

sethaxen added 3 commits January 17, 2023 01:36

Add mcse_sbm to docs

60d6441

Skip high autocorrelation tests for mcse_sbm

2441bcb

Note underestimation for SBM

8e3b06a

sethaxen commented Jan 17, 2023

View reviewed changes

Merge branch 'main' into mcseupdate

e7ca85a

sethaxen added 10 commits January 18, 2023 12:44

Document kwargs passed to mcse

4d55716

Cross-link mcse and ess_rhat docstrings

d9f6734

Document derivation of mcse for std

1c48266

Test type-inferrability of ess_rhat

f072b9e

Make sure ess_rhat for quantiles not promoted

bb47887

Make sure ess_rhat for median type-inferrable

7f61907

Implement specific method for median

a03cc2a

Return missing if any are missing

bac8a3c

Add mcse tests

652b86f

Decrease the number of checks

8dbae84

sethaxen commented Jan 18, 2023

View reviewed changes

sethaxen added 6 commits January 18, 2023 16:32

Make ESS/MCSE for median with with Union{Missing,Real}

d9aff61

Make _fold_around_median type-inferrable

cced4be

Increase tolerance for exhaustive tests

ce9d427

Fix _fold_around_median

787a05f

Fix count of checks

34f3771

Increase the number of draws

d10740a

improves the quality of the estimates and reduces random failures

sethaxen requested a review from devmotion January 18, 2023 20:24

devmotion reviewed Jan 18, 2023

View reviewed changes

sethaxen and others added 6 commits January 18, 2023 22:21

Apply suggestions from code review

cf09e4e

Co-authored-by: David Widmann <[email protected]>

Make sure heideldiag and gewekediag preserve input type

39d74a9

Consistently use first and last for ess_rhat

a575fc8

Copy comment to _fold_around_median

b4eea8d

Make mcse_sbm an internal function

3738347

Update tests

69fe0a1

sethaxen requested a review from devmotion January 18, 2023 22:49

devmotion approved these changes Jan 19, 2023

View reviewed changes

sethaxen merged commit 9e78f78 into main Jan 19, 2023

sethaxen deleted the mcseupdate branch January 19, 2023 15:32

devmotion mentioned this pull request Jan 25, 2023

Stop using Base.@irrational JuliaMath/IrrationalConstants.jl#18

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Redesign of MCSE #63

Redesign of MCSE #63

sethaxen commented Jan 16, 2023

coveralls commented Jan 17, 2023 •

edited

Loading

codecov bot commented Jan 17, 2023 •

edited

Loading

sethaxen Jan 17, 2023

sethaxen Jan 17, 2023 •

edited

Loading

sethaxen Jan 18, 2023

sethaxen Jan 17, 2023

devmotion Jan 18, 2023

sethaxen Jan 18, 2023

sethaxen Jan 18, 2023

sethaxen Jan 18, 2023 •

edited

Loading

devmotion Jan 18, 2023

		SBM tends to underestimate the MCSE, especially for highly autocorrelated chains.
		SBM should only be used as a fallbeck when a specific [`mcse`](@ref) method for

Redesign of MCSE #63

Redesign of MCSE #63

Conversation

sethaxen commented Jan 16, 2023

coveralls commented Jan 17, 2023 • edited Loading

Pull Request Test Coverage Report for Build 3939266728

Warning: This coverage report may be inaccurate.

Details

💛 - Coveralls

codecov bot commented Jan 17, 2023 • edited Loading

Codecov Report

sethaxen Jan 17, 2023

Choose a reason for hiding this comment

sethaxen Jan 17, 2023 • edited Loading

Choose a reason for hiding this comment

sethaxen Jan 18, 2023

Choose a reason for hiding this comment

sethaxen Jan 17, 2023

Choose a reason for hiding this comment

devmotion Jan 18, 2023

Choose a reason for hiding this comment

sethaxen Jan 18, 2023

Choose a reason for hiding this comment

sethaxen Jan 18, 2023

Choose a reason for hiding this comment

sethaxen Jan 18, 2023 • edited Loading

Choose a reason for hiding this comment

devmotion Jan 18, 2023

Choose a reason for hiding this comment

coveralls commented Jan 17, 2023 •

edited

Loading

codecov bot commented Jan 17, 2023 •

edited

Loading

sethaxen Jan 17, 2023 •

edited

Loading

sethaxen Jan 18, 2023 •

edited

Loading