From 8bcdfab589f5748692113bd1ff8ce0f3c20ce74a Mon Sep 17 00:00:00 2001 From: Ray Kim Date: Mon, 20 Nov 2023 20:07:30 +0000 Subject: [PATCH] fix implement rank plots using `RecipesBase` --- Project.toml | 9 ++---- docs/make.jl | 8 ++++- docs/src/example.md | 4 +-- docs/src/general.md | 4 +-- docs/src/index.md | 1 + docs/src/ranksim.md | 10 +++---- src/MCMCTesting.jl | 19 +++--------- ext/MCMCTestingPlotsExt.jl => src/rankplot.jl | 30 +++++++++---------- test/runtests.jl | 4 +-- 9 files changed, 39 insertions(+), 50 deletions(-) rename ext/MCMCTestingPlotsExt.jl => src/rankplot.jl (82%) diff --git a/Project.toml b/Project.toml index b532278..84c38cc 100644 --- a/Project.toml +++ b/Project.toml @@ -7,26 +7,21 @@ HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" MultipleTesting = "f8716d33-7c4a-5097-896f-ce0ecbd3ef6b" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Requires = "ae029012-a4dd-5104-9daa-d747884805df" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -[weakdeps] -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" - -[extensions] -MCMCTestingPlotsExt = "Plots" - [compat] HypothesisTests = "0.11" MultipleTesting = "0.6" ProgressMeter = "1.9" Random = "1" +RecipesBase = "1.3" Requires = "1" StatsBase = "0.34" julia = "1.6" [extras] -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] diff --git a/docs/make.jl b/docs/make.jl index 77d7fc2..016bb51 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,10 +1,16 @@ using MCMCTesting using Documenter +using Plots DocMeta.setdocmeta!(MCMCTesting, :DocTestSetup, :(using MCMCTesting); recursive=true) makedocs(; - modules=[MCMCTesting], + modules=[ + MCMCTesting, + isdefined(Base, :get_extension) ? + Base.get_extension(MCMCTesting, :MCMCTestingPlotsExt) : + RecipesBase.apply_recipe + ], repo="https://github.com/Red-Portal/MCMCTesting.jl/blob/{commit}{path}#{line}", sitename="MCMCTesting.jl", format=Documenter.HTML(; diff --git a/docs/src/example.md b/docs/src/example.md index f585a8c..2ca1000 100644 --- a/docs/src/example.md +++ b/docs/src/example.md @@ -194,8 +194,8 @@ rank_correct = simulate_ranks(test, TestSubject(model, Gibbs()); show_progress=f rank_wrong = simulate_ranks(test, TestSubject(model, GibbsWrongVar()); show_progress=false) param_names = ["θ1 mean", "θ2 mean", "θ1 var", "θ2 var"] -plot(rank_wrong, test; param_names) -plot!(rank_correct, test; param_names) +rankplot(test, rank_wrong; param_names) +rankplot!(test, rank_correct; param_names) savefig("rankplot.svg") nothing ``` diff --git a/docs/src/general.md b/docs/src/general.md index d2919f8..91308db 100644 --- a/docs/src/general.md +++ b/docs/src/general.md @@ -28,13 +28,13 @@ The `model` and `kernel` are then passed to `MCMCTesting` through the following TestSubject ``` -## Simulating a P-Value Through `mcmctest` +## Simulating a P-Value with `mcmctest` Each of the test internally run simulations and compute a single p-value through the following routine: ```@docs mcmctest ``` -## Increasing Power Through `seqmcmctest` +## Increasing Power with `seqmcmctest` `seqmcmctest` (Algorithm 3[^gandyandscott2021]) sequentially calls `mcmctest` to increase the power and ensure a low false rejection rate. Furthermore, the p-values from each component of the statistics are combined through multiple hypothesis adjustment. diff --git a/docs/src/index.md b/docs/src/index.md index d172274..1a90429 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -12,6 +12,7 @@ Unlike simulation-based calibration (SBC; [^TBSVG2018][^YD2023][^MMKBPBHFGV2022] This is because the tests focus on maximizing the power for verify the validity of *individual Markov transitions* instead of a *set of samples*. Furthermore, unlike SBC, the approach of Gandy & Scott[^GS2021] is able to exactly satisfy the assumptions required for the theoretical guarantees. + ## References [^GS2021]: Gandy, A., & Scott, J. (2020). Unit testing for MCMC and other Monte Carlo methods. arXiv preprint arXiv:2001.06465. [^G2004]: Geweke, J. (2004). Getting it right: Joint distribution tests of posterior simulators. Journal of the American Statistical Association, 99(467), 799-804. diff --git a/docs/src/ranksim.md b/docs/src/ranksim.md index 1f2d9ee..1fc4707 100644 --- a/docs/src/ranksim.md +++ b/docs/src/ranksim.md @@ -10,10 +10,11 @@ For more information, refer to the documentation for the [exact rank test](@ref simulate_ranks ``` -## Visualizing Ranks +## Visualizing Ranks with `Plots` We provide a `Plots` recipe for visualizing the ranks: -```@docs -MCMCTesting.plotranks + +```docs +rankplot ``` This can be used as follows: @@ -24,10 +25,9 @@ using MCMCTesting # Set up the simulation ranks = simulate_ranks(test, subject) -plot(ranks, test; param_names) +rankplot(test, ranks; param_names) ``` Also refer to the [tutorial](@ref tutorial) for a working example. - ## References [^GS2021]: Gandy, A., & Scott, J. (2020). Unit testing for MCMC and other Monte Carlo methods. arXiv preprint arXiv:2001.06465. diff --git a/src/MCMCTesting.jl b/src/MCMCTesting.jl index db01a24..c212fdc 100644 --- a/src/MCMCTesting.jl +++ b/src/MCMCTesting.jl @@ -11,15 +11,16 @@ export markovchain_transition, mcmctest, seqmcmctest, - simulate_ranks + simulate_ranks, + RankPlot using HypothesisTests using MultipleTesting using ProgressMeter using Random +using RecipesBase using StatsBase - """ sample_joint(rng, model) @@ -120,10 +121,8 @@ struct TestSubject{M, K} kernel::K end - abstract type AbstractMCMCTest end - function markovchain_multiple_transition( rng::Random.AbstractRNG, model, kernel, n_steps::Int, θ, y ) @@ -137,16 +136,6 @@ include("defaults.jl") include("twosampletest.jl") include("exactranktest.jl") include("seqtest.jl") +include("rankplot.jl") -if !isdefined(Base, :get_extension) - using Requires -end - -@static if !isdefined(Base, :get_extension) - function __init__() - @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin - include("../ext/MCMCTestingPlotsExt.jl") - end - end -end end diff --git a/ext/MCMCTestingPlotsExt.jl b/src/rankplot.jl similarity index 82% rename from ext/MCMCTestingPlotsExt.jl rename to src/rankplot.jl index 679018b..a4000f8 100644 --- a/ext/MCMCTestingPlotsExt.jl +++ b/src/rankplot.jl @@ -1,28 +1,27 @@ -module MCMCTestingPlotsExt - -if isdefined(Base, :get_extension) - using Plots - using MCMCTesting -else - using ..Plots - using ..MCMCTesting -end - """ - plotranks(ranks, test; kwargs...) + rankplot(test, ranks; kwargs...) -`Plots` recipe for plotting the simulated ranks from `simulate_ranks`. -`Plots` must be imported *before* `MCMCTesting` to use this plot recipe. +Plot the simulated ranks using `simulate_ranks`. + +!!! info + `Plots` must be imported *before* `MCMCTesting` to use this plot recipe. # Arguments -- ranks: The output of `simulate_rank`. - test::ExactRankTest: The exact rank test object used to simulate the ranks. +- ranks: The output of `simulate_rank`. # Keyword Arguments - stats_names: The name for the statistics used in the rank simulation. The default argument automatically assign default names. (Default: :auto). +- Keyword arguments corresponding `Plots` attributes, such as `bins`, `layout`, `size`, may apply. """ -@recipe function plotranks(ranks, test::ExactRankTest; stat_names=:auto) +@userplot RankPlot +@recipe function f(h::RankPlot; stat_names = :auto) + if length(h.args) != 2 || !(typeof(h.args[1]) <: ExactRankTest) + error("rankplot should be given a `<: ExctRankTest` as first argument. Got: $(typeof(h.args)).") + end + test, ranks = h.args + n_max_rank = test.n_mcmc_steps n_samples = test.n_samples binprob = 1/n_max_rank @@ -93,4 +92,3 @@ end end end end -end diff --git a/test/runtests.jl b/test/runtests.jl index 8d0ec7b..4c0a49b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -89,11 +89,11 @@ include("models/normalnormalgibbs.jl") ranks = simulate_ranks(test, subject; show_progress=false) @testset "default arguments" begin - plot(ranks, test) + rankplot(test, ranks) end @testset "custom param names" begin - plot(ranks, test; param_names=["test$(idx)" for idx in size(ranks, 1)]) + rankplot(test, ranks; param_names=["test$(idx)" for idx in size(ranks, 1)]) end end end