Skip to content

Commit

Permalink
fix implement rank plots using RecipesBase
Browse files Browse the repository at this point in the history
  • Loading branch information
Red-Portal committed Nov 20, 2023
1 parent 8392faa commit 8bcdfab
Show file tree
Hide file tree
Showing 9 changed files with 39 additions and 50 deletions.
9 changes: 2 additions & 7 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
8 changes: 7 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
@@ -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(;
Expand Down
4 changes: 2 additions & 2 deletions docs/src/example.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
```
Expand Down
4 changes: 2 additions & 2 deletions docs/src/general.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down
1 change: 1 addition & 0 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
10 changes: 5 additions & 5 deletions docs/src/ranksim.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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.
19 changes: 4 additions & 15 deletions src/MCMCTesting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
)
Expand All @@ -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
30 changes: 14 additions & 16 deletions ext/MCMCTestingPlotsExt.jl → src/rankplot.jl
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -93,4 +92,3 @@ end
end
end
end
end
4 changes: 2 additions & 2 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 8bcdfab

Please sign in to comment.