Extract resampling methods

Project.toml

@@ -3,5 +3,9 @@ uuid = "576499cb-2369-40b2-a588-c64705576edc"
 authors = ["TuringLang"]
 version = "0.1.0"
 
+[deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+
 [compat]
+Distributions = "0.23, 0.24"
 julia = "1.3"

src/AdvancedPS.jl

@@ -1,5 +1,6 @@
 module AdvancedPS
 
-# Write your package code here.
+import Distributions
 
+include("resampling.jl")
 end

src/container.jl

@@ -334,15 +334,6 @@ function sweep!(pc::ParticleContainer, resampler)
  return logevidence
 end
 
-struct ResampleWithESSThreshold{R, T<:Real}
- resampler::R
- threshold::T
-end
-
-function ResampleWithESSThreshold(resampler = Turing.Inference.resample_systematic)
- ResampleWithESSThreshold(resampler, 0.5)
-end
-
 function resample_propagate!(
  pc::ParticleContainer,
  resampler::ResampleWithESSThreshold,

src/resampling.jl

@@ -0,0 +1,159 @@
+####
+#### Resampling schemes for particle filters
+####
+
+# Some references
+# - http://arxiv.org/pdf/1301.4019.pdf
+# - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
+# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering
+
+struct ResampleWithESSThreshold{R, T<:Real}
+ resampler::R
+ threshold::T
+end
+
+function ResampleWithESSThreshold(resampler = resample)
+ ResampleWithESSThreshold(resampler, 0.5)
+end
+
+# Default resampling scheme
+function resample(w::AbstractVector{<:Real}, num_particles::Integer=length(w))
+ return resample_systematic(w, num_particles)
+end
+
+# More stable, faster version of rand(Categorical)
+function randcat(p::AbstractVector{<:Real})
+ T = eltype(p)
+ r = rand(T)
+ s = 1
+ for j in eachindex(p)
+ r -= p[j]
+ if r <= zero(T)
+ s = j
+ break
+ end
+ end
+ return s
+end
+
+function resample_multinomial(w::AbstractVector{<:Real}, num_particles::Integer)
+ return rand(Distributions.sampler(Distributions.Categorical(w)), num_particles)
+end
+
+function resample_residual(w::AbstractVector{<:Real}, num_particles::Integer)
+ # Pre-allocate array for resampled particles
+ indices = Vector{Int}(undef, num_particles)
+
+ # deterministic assignment
+ residuals = similar(w)
+ i = 1
+ @inbounds for j in 1:length(w)
+ x = num_particles * w[j]
+ floor_x = floor(Int, x)
+ for k in 1:floor_x
+ indices[i] = j
+ i += 1
+ end
+ residuals[j] = x - floor_x
+ end
+
+ # sampling from residuals
+ if i <= num_particles
+ residuals ./= sum(residuals)
+ rand!(Distributions.Categorical(residuals), view(indices, i:num_particles))
+ end
+
+ return indices
+end
+
+
+"""
+ resample_stratified(weights, n)
+
+Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
+generated by stratified resampling.
+
+In stratified resampling `n` ordered random numbers `u₁`, ..., `uₙ` are generated, where
+``uₖ \\sim U[(k - 1) / n, k / n)``. Based on these numbers the samples `x₁`, ..., `xₙ`
+are selected according to the multinomial distribution defined by the normalized `weights`,
+i.e., `xᵢ = j` if and only if
+``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
+"""
+function resample_stratified(weights::AbstractVector{<:Real}, n::Integer)
+ # check input
+ m = length(weights)
+ m > 0 || error("weight vector is empty")
+
+ # pre-calculations
+ @inbounds v = n * weights[1]
+
+ # generate all samples
+ samples = Array{Int}(undef, n)
+ sample = 1
+ @inbounds for i in 1:n
+ # sample next `u` (scaled by `n`)
+ u = oftype(v, i - 1 + rand())
+
+ # as long as we have not found the next sample
+ while v < u
+ # increase and check the sample
+ sample += 1
+ sample > m &&
+ error("sample could not be selected (are the weights normalized?)")
+
+ # update the cumulative sum of weights (scaled by `n`)
+ v += n * weights[sample]
+ end
+
+ # save the next sample
+ samples[i] = sample
+ end
+
+ return samples
+end
+
+"""
+ resample_systematic(weights, n)
+
+Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
+generated by systematic resampling.
+
+In systematic resampling a random number ``u \\sim U[0, 1)`` is used to generate `n` ordered
+numbers `u₁`, ..., `uₙ` where ``uₖ = (u + k − 1) / n``. Based on these numbers the samples
+`x₁`, ..., `xₙ` are selected according to the multinomial distribution defined by the
+normalized `weights`, i.e., `xᵢ = j` if and only if
+``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
+"""
+function resample_systematic(weights::AbstractVector{<:Real}, n::Integer)
+ # check input
+ m = length(weights)
+ m > 0 || error("weight vector is empty")
+
+ # pre-calculations
+ @inbounds v = n * weights[1]
+ u = oftype(v, rand())
+
+ # find all samples
+ samples = Array{Int}(undef, n)
+ sample = 1
+ @inbounds for i in 1:n
+ # as long as we have not found the next sample
+ while v < u
+ # increase and check the sample
+ sample += 1
+ sample > m &&
+ error("sample could not be selected (are the weights normalized?)")
+
+ # update the cumulative sum of weights (scaled by `n`)
+ v += n * weights[sample]
+ end
+
+ # save the next sample
+ samples[i] = sample
+
+ # update `u`
+ u += one(u)
+ end
+
+ return samples
+end

src/smc.jl

@@ -348,155 +348,4 @@ end
 function DynamicPPL.observe(spl::Sampler{<:Union{PG,SMC}}, dist::Distribution, value, vi)
  produce(logpdf(dist, value))
  return 0
-end
-
-####
-#### Resampling schemes for particle filters
-####
-
-# Some references
-# - http://arxiv.org/pdf/1301.4019.pdf
-# - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
-# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering
-
-# Default resampling scheme
-function resample(w::AbstractVector{<:Real}, num_particles::Integer=length(w))
- return resample_systematic(w, num_particles)
-end
-
-# More stable, faster version of rand(Categorical)
-function randcat(p::AbstractVector{<:Real})
- T = eltype(p)
- r = rand(T)
- s = 1
- for j in eachindex(p)
- r -= p[j]
- if r <= zero(T)
- s = j
- break
- end
- end
- return s
-end
-
-function resample_multinomial(w::AbstractVector{<:Real}, num_particles::Integer)
- return rand(Distributions.sampler(Categorical(w)), num_particles)
-end
-
-function resample_residual(w::AbstractVector{<:Real}, num_particles::Integer)
- # Pre-allocate array for resampled particles
- indices = Vector{Int}(undef, num_particles)
-
- # deterministic assignment
- residuals = similar(w)
- i = 1
- @inbounds for j in 1:length(w)
- x = num_particles * w[j]
- floor_x = floor(Int, x)
- for k in 1:floor_x
- indices[i] = j
- i += 1
- end
- residuals[j] = x - floor_x
- end
-
- # sampling from residuals
- if i <= num_particles
- residuals ./= sum(residuals)
- rand!(Categorical(residuals), view(indices, i:num_particles))
- end
-
- return indices
-end
-
-
-"""
- resample_stratified(weights, n)
-
-Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
-generated by stratified resampling.
-
-In stratified resampling `n` ordered random numbers `u₁`, ..., `uₙ` are generated, where
-``uₖ \\sim U[(k - 1) / n, k / n)``. Based on these numbers the samples `x₁`, ..., `xₙ`
-are selected according to the multinomial distribution defined by the normalized `weights`,
-i.e., `xᵢ = j` if and only if
-``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
-"""
-function resample_stratified(weights::AbstractVector{<:Real}, n::Integer)
- # check input
- m = length(weights)
- m > 0 || error("weight vector is empty")
-
- # pre-calculations
- @inbounds v = n * weights[1]
-
- # generate all samples
- samples = Array{Int}(undef, n)
- sample = 1
- @inbounds for i in 1:n
- # sample next `u` (scaled by `n`)
- u = oftype(v, i - 1 + rand())
-
- # as long as we have not found the next sample
- while v < u
- # increase and check the sample
- sample += 1
- sample > m &&
- error("sample could not be selected (are the weights normalized?)")
-
- # update the cumulative sum of weights (scaled by `n`)
- v += n * weights[sample]
- end
-
- # save the next sample
- samples[i] = sample
- end
-
- return samples
-end
-
-"""
- resample_systematic(weights, n)
-
-Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
-generated by systematic resampling.
-
-In systematic resampling a random number ``u \\sim U[0, 1)`` is used to generate `n` ordered
-numbers `u₁`, ..., `uₙ` where ``uₖ = (u + k − 1) / n``. Based on these numbers the samples
-`x₁`, ..., `xₙ` are selected according to the multinomial distribution defined by the
-normalized `weights`, i.e., `xᵢ = j` if and only if
-``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
-"""
-function resample_systematic(weights::AbstractVector{<:Real}, n::Integer)
- # check input
- m = length(weights)
- m > 0 || error("weight vector is empty")
-
- # pre-calculations
- @inbounds v = n * weights[1]
- u = oftype(v, rand())
-
- # find all samples
- samples = Array{Int}(undef, n)
- sample = 1
- @inbounds for i in 1:n
- # as long as we have not found the next sample
- while v < u
- # increase and check the sample
- sample += 1
- sample > m &&
- error("sample could not be selected (are the weights normalized?)")
-
- # update the cumulative sum of weights (scaled by `n`)
- v += n * weights[sample]
- end
-
- # save the next sample
- samples[i] = sample
-
- # update `u`
- u += one(u)
- end
-
- return samples
-end
+end

test/resampling.jl

@@ -0,0 +1,17 @@
+@testset "resampling.jl" begin
+ D = [0.3, 0.4, 0.3]
+ num_samples = Int(1e6)
+
+ resSystematic = AdvancedPS.resample_systematic(D, num_samples )
+ resStratified = AdvancedPS.resample_stratified(D, num_samples )
+ resMultinomial= AdvancedPS.resample_multinomial(D, num_samples )
+ resResidual = AdvancedPS.resample_residual(D, num_samples )
+ AdvancedPS.resample(D)
+ resSystematic2= AdvancedPS.resample(D, num_samples )
+
+ @test sum(resSystematic .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples
+ @test sum(resSystematic2 .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples
+ @test sum(resStratified .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples
+ @test sum(resMultinomial .== 2) ≈ (num_samples * 0.4) atol=1e-2*num_samples
+ @test sum(resResidual .== 2) ≈ (num_samples * 0.4) atol=1e-2*num_samples
+end

test/runtests.jl

@@ -2,5 +2,5 @@ using AdvancedPS
 using Test
 
 @testset "AdvancedPS.jl" begin
- # Write your tests here.
+ @testset "Resampling tests" begin include("resampling.jl") end
 end

test/smc.jl

@@ -235,19 +235,3 @@ end
 # end
 # end
 
-@turing_testset "resample.jl" begin
- D = [0.3, 0.4, 0.3]
- num_samples = Int(1e6)
- resSystematic = Turing.Inference.resample_systematic(D, num_samples )
- resStratified = Turing.Inference.resample_stratified(D, num_samples )
- resMultinomial= Turing.Inference.resample_multinomial(D, num_samples )
- resResidual = Turing.Inference.resample_residual(D, num_samples )
- Turing.Inference.resample(D)
- resSystematic2=Turing.Inference.resample(D, num_samples )
-
- @test sum(resSystematic .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples
- @test sum(resSystematic2 .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples
- @test sum(resStratified .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples
- @test sum(resMultinomial .== 2) ≈ (num_samples * 0.4) atol=1e-2*num_samples
- @test sum(resResidual .== 2) ≈ (num_samples * 0.4) atol=1e-2*num_samples
-end