From e674961e12b9a6d3d4c38ee32ef537e69944f776 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 2 Dec 2020 13:20:16 +0100 Subject: [PATCH] Extract resampling methods --- Project.toml | 4 ++ src/AdvancedPS.jl | 3 +- src/container.jl | 9 --- src/resampling.jl | 159 +++++++++++++++++++++++++++++++++++++++++++++ src/smc.jl | 153 +------------------------------------------ test/resampling.jl | 17 +++++ test/runtests.jl | 2 +- test/smc.jl | 16 ----- 8 files changed, 184 insertions(+), 179 deletions(-) create mode 100644 src/resampling.jl create mode 100644 test/resampling.jl diff --git a/Project.toml b/Project.toml index e10e01bc..adca95de 100644 --- a/Project.toml +++ b/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" diff --git a/src/AdvancedPS.jl b/src/AdvancedPS.jl index a86133e1..8fcbba3d 100644 --- a/src/AdvancedPS.jl +++ b/src/AdvancedPS.jl @@ -1,5 +1,6 @@ module AdvancedPS -# Write your package code here. +import Distributions +include("resampling.jl") end diff --git a/src/container.jl b/src/container.jl index 35353f01..4f4c9822 100644 --- a/src/container.jl +++ b/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, diff --git a/src/resampling.jl b/src/resampling.jl new file mode 100644 index 00000000..1cb0ef40 --- /dev/null +++ b/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 diff --git a/src/smc.jl b/src/smc.jl index f2354b4d..ac5a5c11 100644 --- a/src/smc.jl +++ b/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 \ No newline at end of file diff --git a/test/resampling.jl b/test/resampling.jl new file mode 100644 index 00000000..ee1ff953 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index d801d119..067506a3 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/smc.jl b/test/smc.jl index a987cc81..ec2bca09 100644 --- a/test/smc.jl +++ b/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