Basic Filtering Structure

examples/particle-mcmc/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001"
+AdvancedMH = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+GaussianDistributions = "43dcc890-d446-5863-8d1a-14597580bb8d"
+PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+SSMProblems = "26aad666-b158-4e64-9d35-0e672562fa48"
+StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

examples/particle-mcmc/script.jl

@@ -0,0 +1,85 @@
+using AdvancedMH
+using CairoMakie
+
+include("simple-filters.jl")
+
+true_params, simulation_model = let T = Float32
+    θ = randexp(T, 3)
+    dyn = LinearGaussianLatentDynamics(T[1 1;0 1], diagm(θ[1:2]))
+    obs = LinearGaussianObservationProcess(T[0.5 0.5], diagm(θ[3:end]))
+    θ, StateSpaceModel(dyn, obs)
+end
+
+# simulate data
+rng = MersenneTwister(1234)
+_, _, data = sample(rng, simulation_model, 150)
+
+# consider a default Gamma prior with Float32s
+prior_dist = product_distribution(Gamma(1f0), Gamma(1f0), Gamma(1f0))
+
+# test the adaptive resampling procedure
+sample(rng, simulation_model, data, BF(512, 0.1); debug=true);
+
+
+#=
+    Not crazy about this structure, especially since the RNG is referenced on
+    the global scope. I think we can make a PMCMC sampler type which includes
+    the filter algorithm within the sampler definition.
+
+    Another issue is that we lose information on the states. Granted, this is
+    also by design since that would cost a considerable amount of memory, but
+    is useful nonetheless. This also needs to interface with bundle_samples()
+    different than ususal, since we have the parameter space and the filtered
+    states.
+=#
+function density(θ::Vector{T}) where {T<:Real}
+    if insupport(prior_dist, θ)
+        dyn = LinearGaussianLatentDynamics(T[1 1;0 1], diagm(θ[1:2]))
+        obs = LinearGaussianObservationProcess(T[0.5 0.5], diagm(θ[3:end]))
+
+        # _, ll = sample(rng, StateSpaceModel(dyn, obs), data, BF(512))
+        _, ll = sample(rng, StateSpaceModel(dyn, obs), data, KF())
+        return ll + logpdf(prior_dist, θ)
+    else
+        return -Inf
+    end
+end
+
+# plug it into the DensityModel interface for now
+pmmh  = RWMH(MvNormal(zeros(Float32, 3), (0.01f0)*I))
+model = DensityModel(density)
+
+# works with AdvancedMH out of the box
+chains = sample(model, pmmh, 10_000)
+burn_in = 1_000
+
+# plot the posteriors
+hist_plots = begin
+    param_post = hcat(getproperty.(chains[burn_in:end], :params)...)
+    fig = Figure(size = (1200, 400))
+
+    for i in 1:3
+        # plot the posteriors with burn-in
+        hist(
+            fig[1, i],
+            param_post[i, :],
+            color = :gray,
+            strokewidth = 1,
+            normalization = :pdf
+        )
+
+        # plot the true values
+        vlines!(
+            fig[1, i],
+            true_params[i],
+            color = :red,
+            linestyle = :dash,
+            linewidth = 2
+        )
+    end
+
+    fig
+end
+
+# this is useful for SMC algorithms like SMC² or density tempered SMC
+acc_ratio = mean(getproperty.(chains, :accepted))