Merge pull request #50 from luiarthur/master

Implemented `thinning` option in `mcmcsample`

Project.toml

@@ -3,7 +3,7 @@ uuid = "80f14c24-f653-4e6a-9b94-39d6b0f70001"
 keywords = ["markov chain monte carlo", "probablistic programming"]
 license = "MIT"
 desc = "A lightweight interface for common MCMC methods."
-version = "2.1.0"
+version = "2.2.0"
 
 [deps]
 BangBang = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"

README.md

@@ -72,6 +72,7 @@ are:
   `callback(rng, model, sampler, sample, iteration)` is called after every sampling step,
   where `sample` is the most recent sample of the Markov chain and `iteration` is the current iteration
 - `discard_initial` (default: `0`): number of initial samples that are discarded
+- `thinning` (default: `1`): factor by which to thin samples.
 
 Additionally, AbstractMCMC defines the abstract type `AbstractChains` for Markov chains and the
 method `AbstractMCMC.chainscat(::AbstractChains...)` for concatenating multiple chains.

src/sample.jl

@@ -65,12 +65,13 @@ function mcmcsample(
     progressname = "Sampling",
     callback = nothing,
     discard_initial = 0,
+    thinning = 1,
     chain_type::Type=Any,
     kwargs...
 )
     # Check the number of requested samples.
     N > 0 || error("the number of samples must be ≥ 1")
-    Ntotal = N + discard_initial
+    Ntotal = thinning * (N - 1) + discard_initial + 1
 
     @ifwithprogresslogger progress name=progressname begin
         # Obtain the initial sample and state.
@@ -96,7 +97,20 @@ function mcmcsample(
         progress && ProgressLogging.@logprogress (1 + discard_initial) / Ntotal
 
         # Step through the sampler.
+        itotal = 1 + discard_initial
         for i in 2:N
+            # Discard thinned samples.
+            for _ in 1:(thinning - 1)
+                # Obtain the next sample and state.
+                sample, state = step(rng, model, sampler, state; kwargs...)
+
+                # Update progress bar.
+                if progress
+                    itotal += 1
+                    ProgressLogging.@logprogress itotal / Ntotal
+                end
+            end
+
             # Obtain the next sample and state.
             sample, state = step(rng, model, sampler, state; kwargs...)
 
@@ -107,7 +121,10 @@ function mcmcsample(
             samples = save!!(samples, sample, i, model, sampler, N; kwargs...)
 
             # Update the progress bar.
-            progress && ProgressLogging.@logprogress (i + discard_initial) / Ntotal
+            if progress
+                itotal += 1
+                ProgressLogging.@logprogress itotal / Ntotal
+            end
         end
     end
 
@@ -141,6 +158,7 @@ function mcmcsample(
     progressname = "Convergence sampling",
     callback = nothing,
     discard_initial = 0,
+    thinning = 1,
     kwargs...
 )
     @ifwithprogresslogger progress name=progressname begin
@@ -164,6 +182,12 @@ function mcmcsample(
         i = 2
 
         while !isdone(rng, model, sampler, samples, i; progress=progress, kwargs...)
+            # Discard thinned samples.
+            for _ in 1:(thinning - 1)
+                # Obtain the next sample and state.
+                sample, state = step(rng, model, sampler, state; kwargs...)
+            end
+
             # Obtain the next sample and state.
             sample, state = step(rng, model, sampler, state; kwargs...)
 

test/sample.jl

@@ -226,6 +226,22 @@
         @test !ismissing(chain[1].a)
     end
 
+    @testset "Thin chain by a factor of `thinning`" begin
+        # Run a thinned chain with `N` samples thinned by factor of `thinning`.
+        Random.seed!(1234)
+        N = 100
+        thinning = 3
+        chain = sample(MyModel(), MySampler(), N; sleepy = true, thinning = thinning)
+        @test length(chain) == N
+        @test ismissing(chain[1].a)
+
+        # Repeat sampling without thinning.
+        Random.seed!(1234)
+        ref_chain = sample(MyModel(), MySampler(), N * thinning; sleepy = true)
+        @test all(chain[i].a === ref_chain[(i - 1) * thinning + 1].a for i in 1:N)
+    end
+
+
     @testset "Sample without predetermined N" begin
         Random.seed!(1234)
         chain = sample(MyModel(), MySampler())
@@ -238,6 +254,12 @@
         bmean = mean(x.b for x in chain)
         @test !ismissing(chain[1].a)
         @test abs(bmean) <= 0.001 && length(chain) < 10_000
+
+        # Thin chain by a factor of `thinning`.
+        chain = sample(MyModel(), MySampler(); thinning = 3)
+        bmean = mean(x.b for x in chain)
+        @test ismissing(chain[1].a)
+        @test abs(bmean) <= 0.001 && length(chain) < 10_000
     end
 
     @testset "Sample vector of `NamedTuple`s" begin