update documentation with correct Gibbs sampler

docs/Project.toml

@@ -3,6 +3,7 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 MCMCTesting = "9963b6a1-5d46-439c-8efc-3a487843c7fa"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]

docs/src/exactranktest.md

@@ -2,7 +2,7 @@
 # [Exact Rank Hypothesis Tests](@id exactrank)
 
 ## Introduction
-The exact rank hypothesis testing strategy is based on Algorithm 2 by Gandy & Scott (2021)[^gs2021].
+The exact rank hypothesis testing strategy is based on Algorithm 2 by Gandy & Scott (2021)[^GS2021].
 
 ## `ExactRankTest`
 The ranks are computed by simulating a single Markov chain backwards and forward.

docs/src/example.md

@@ -2,11 +2,12 @@
 # [Tutorial](@id tutorial)
 
 ## Problem Setup
-Let's consider a simple Normal-Normal model with a shared scale:
+Let's consider a simple Normal-Normal model:
 ```math
 \begin{aligned}
-\theta &\sim \mathrm{normal}\left(0,   \sigma^2\right) \\
-y      &\sim \mathrm{normal}\left(\mu, \sigma^2\right)
+\theta_1 &\sim \mathrm{normal}\left(0,   \sigma^2\right) \\
+\theta_2 &\sim \mathrm{normal}\left(0,   \sigma^2\right) \\
+y        &\sim \mathrm{normal}\left(\theta_1 + \theta_2, \sigma_{\epsilon}^2\right).
 \end{aligned}
 ```
 
@@ -16,44 +17,30 @@ using Random
 using Distributions
 
 struct Model
-    σ::Float64
-    y::Float64
-end
-
-function logdensity(model::Model, θ)
-    σ, y = model.σ, model.y
-    logpdf(Normal(0, σ), only(θ)) + logpdf(Normal(0, σ), y)
+    sigma    ::Float64
+    sigma_eps::Float64
+    y        ::Float64
 end
 nothing
 ```
 For sampling from the posterior, a simple Gibbs sampling strategy is possible:
 ```@example started
 struct Gibbs end
 
-function step(rng::Random.AbstractRNG, model::Model, ::Gibbs, θ)
-    y = model.y
-    σ = model.σ
-    rand(rng, MvNormal([y]/2, σ))
-end
-nothing
-```
-We could also use the classic symmetric random walk Metropolis-Hastings sampler:
-```@example started
-struct RWMH
-    σ::Float64
+function complete_conditional(θ::Real, σ²::Real, σ²_ϵ::Real, y::Real)
+    μ = σ²/(σ²_ϵ + σ²)*(y - θ)
+    σ = 1/sqrt(1/σ²_ϵ + 1/σ²)
+    Normal(μ, σ)
 end
 
-function step(rng::Random.AbstractRNG, model::Model, kernel::RWMH, θ)
-    σ = kernel.σ
-    θ′ = rand(rng, MvNormal(θ, σ))
-    ℓπ = logdensity(model, θ)
-    ℓπ′ = logdensity(model, θ′)
-    ℓα = ℓπ′ - ℓπ
-    if log(rand(rng)) ≤ ℓα
-        θ′
-    else
-        θ
-    end
+function step(rng::Random.AbstractRNG, model::Model, ::Gibbs, θ)
+    θ    = copy(θ)
+    y    = model.y
+    σ²   = model.sigma^2
+    σ²_ϵ = model.sigma_eps^2
+    θ[1] = rand(rng, complete_conditional(θ[2], σ², σ²_ϵ, y))
+    θ[2] = rand(rng, complete_conditional(θ[1], σ², σ²_ϵ, y))
+    θ
 end
 nothing
 ```
@@ -65,16 +52,16 @@ This is done by as follows:
 using MCMCTesting
 
 function MCMCTesting.sample_joint(rng::Random.AbstractRNG, model::Model)
-    σ = model.σ
-    θ = rand(rng, Normal(0, σ))
-    y = rand(rng, Normal(θ, σ))
-    [θ], [y]
+    θ₁ = rand(rng, Normal(0, model.sigma))
+    θ₂ = rand(rng, Normal(0, model.sigma))
+    θ  = [θ₁, θ₂]
+    y  = rand(rng, Normal(θ[1] + θ[2], model.sigma_eps))
+    θ, y
 end
 nothing
 ```
 
-The Gibbs sampler can be connected to `MCMCTesting` by implementing the following:
-
+The Gibbs samplers can be connected to `MCMCTesting` by implementing the following:
 ```@example started
 using Accessors
 
@@ -86,33 +73,59 @@ function MCMCTesting.markovchain_transition(
 end
 nothing
 ```
-Let's check that the implementation is correct by 
 
+Let's check that the implementation is correct by 
 ```@example started
-model = Model(1., 1.)
+model = Model(1., .5, randn())
 kernel = Gibbs()
 test = TwoSampleTest(100, 100)
 subject = TestSubject(model, kernel)
 seqmcmctest(test, subject, 0.0001, 100; show_progress=false)
 ```
 `true` means that the tests have passed.
-Now, let's consider two erroneous implementations:
 
+Now, let's consider two erroneous implementations:
 ```@example started
 struct GibbsWrongMean end
 
+function complete_conditional_wrongmean(θ::Real, σ²::Real, σ²_ϵ::Real, y::Real)
+    μ = σ²/(σ²_ϵ + σ²)*(y + θ)
+    σ = 1/sqrt(1/σ²_ϵ + 1/σ²)
+    Normal(μ, σ)
+end
+
 function step(rng::Random.AbstractRNG, model::Model, ::GibbsWrongMean, θ)
-    y = model.y
-    σ = model.σ
-    rand(rng, MvNormal([y], σ/2))
+    θ    = copy(θ)
+    y    = model.y
+    σ²   = model.sigma^2
+    σ²_ϵ = model.sigma_eps^2
+
+    θ[1] = rand(rng, complete_conditional_wrongmean(θ[2], σ², σ²_ϵ, y))
+    θ[2] = rand(rng, complete_conditional_wrongmean(θ[1], σ², σ²_ϵ, y))
+    θ
 end
 
 struct GibbsWrongVar end
 
+function complete_conditional_wrongvar(θ::Real, σ²::Real, σ²_ϵ::Real, y::Real)
+    μ = σ²/(σ²_ϵ + σ²)*(y - θ)
+    Normal(μ, sqrt(σ²))
+end
+
 function step(rng::Random.AbstractRNG, model::Model, ::GibbsWrongVar, θ)
-    y = model.y
-    σ = model.σ
-    rand(rng, MvNormal([y/2], 2*σ))
+    θ    = copy(θ)
+    y    = model.y
+    σ²   = model.sigma^2
+    σ²_ϵ = model.sigma_eps^2
+
+    if rand(Bernoulli(0.5))
+        θ[1] = rand(rng, complete_conditional_wrongvar(θ[2], σ², σ²_ϵ, y))
+        θ[2] = rand(rng, complete_conditional_wrongvar(θ[1], σ², σ²_ϵ, y))
+    else
+        θ[2] = rand(rng, complete_conditional_wrongvar(θ[1], σ², σ²_ϵ, y))
+        θ[1] = rand(rng, complete_conditional_wrongvar(θ[2], σ², σ²_ϵ, y))
+    end
+    θ
 end
 nothing
 ```
@@ -129,4 +142,62 @@ subject = TestSubject(model, kernel)
 seqmcmctest(test, subject, 0.0001, 100; show_progress=false)
 ```
 
-## Visualizing the ranks
+## Visualizing Simulated Ranks
+`MCMCTesting` also provides some basic plot recipes for visualizing the simulated rank for the [exact rank test](@ref exactrank).
+For this, `Plots` must be imported *before* `MCMCTesting`.
+
+```@example started
+using Plots
+gr()
+using MCMCTesting
+nothing
+```
+
+Also, since the exact rank test explicitly requires the MCMC kernel to be reversible, we modify the previous Gibbs sampler to use a random scan order.
+```@example started
+function step(rng::Random.AbstractRNG, model::Model, ::Gibbs, θ)
+    θ    = copy(θ)
+    y    = model.y
+    σ²   = model.sigma^2
+    σ²_ϵ = model.sigma_eps^2
+    if rand(Bernoulli(0.5))
+        θ[1] = rand(rng, complete_conditional(θ[2], σ², σ²_ϵ, y))
+        θ[2] = rand(rng, complete_conditional(θ[1], σ², σ²_ϵ, y))
+    else
+        θ[2] = rand(rng, complete_conditional(θ[1], σ², σ²_ϵ, y))
+        θ[1] = rand(rng, complete_conditional(θ[2], σ², σ²_ϵ, y))
+    end
+    θ
+end
+
+function step(rng::Random.AbstractRNG, model::Model, ::GibbsWrongVar, θ)
+    θ    = copy(θ)
+    y    = model.y
+    σ²   = model.sigma^2
+    σ²_ϵ = model.sigma_eps^2
+    if rand(Bernoulli(0.5))
+        θ[1] = rand(rng, complete_conditional_wrongvar(θ[2], σ², σ²_ϵ, y))
+        θ[2] = rand(rng, complete_conditional_wrongvar(θ[1], σ², σ²_ϵ, y))
+    else
+        θ[2] = rand(rng, complete_conditional_wrongvar(θ[1], σ², σ²_ϵ, y))
+        θ[1] = rand(rng, complete_conditional_wrongvar(θ[2], σ², σ²_ϵ, y))
+    end
+    θ
+end
+nothing
+```
+Then, we can simulate the ranks and then plot them using `Plots.
+```@example started
+test = ExactRankTest(1000, 30, 10)
+
+rank_correct = simulate_ranks(test, TestSubject(model, Gibbs()); show_progress=false)
+rank_wrong = simulate_ranks(test, TestSubject(model, GibbsWrongVar()); show_progress=false)
+
+plot(rank_wrong[1:1,:],    test, alpha=0.3)
+plot!(rank_correct[1:1,:], test, alpha=0.3)
+savefig("rankplot.svg")
+nothing
+```
+![](rankplot.svg)
+
+We can see that the ranks of the erroneous kernel are not uniform.