Best library for particle filtering from a generative model?

[JasonPekos]

Hi! I’m looking for a library that will let me perform particle filtering for a generative model of a markov process, similar to the incredible pomp in R:

https://kingaa.github.io/pomp/

or particles in Python.

Basically, I want to specify some sort of stochastic generator f such that:

x_t = f(x_{t-1})

and some sort of measurement model for my data y_t:

y_t \sim \operatorname{Binomal}(x_t, p)

and then do something like:

particle_filter(generator, measurement, particles = 1000)

or similar. Maybe Ideally it would work with the SciML syntax for the process model f(x_{t-1}), so something like:

# from julia epi-recipes
function sir_markov!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    ifrac = rate_to_proportion(β*c*I/N,δt)
    rfrac = rate_to_proportion(γ,δt)
    infection=rand(Binomial(S,ifrac))
    recovery=rand(Binomial(I,rfrac))
    @inbounds begin
        du[1] = S-infection
        du[2] = I+infection-recovery
        du[3] = R+recovery
    end
    nothing
end;

Just so the interface is standardized and I don’t run into a giant headache of uninformative stacktraces trying to write a model that the solver will accept.

---

The lay of the land:

I’ve taken a quick look at the packages that seem like good candidates, but I’m not sure what the best option is.

1. ParticleFilters.jl: is a subpackage of POMDPs.jl. This looks like it might be the best candidate, but I’m not sure if there’s a more focused standalone package I’m missing. Most of the documentation here is in the main package, and relates specifically to decision processes, which I’m not familiar with. I get the sense that this is at least a little bit of an internal subpackage for the main POMDPs.jl interface.

2. LowLevelParticleFilters.jl Also looks interesting, but maybe a bit more involved to set up.

3. GenParticleFilters.jl Part of Gen.

4. AdvancedPS.jl Part of Turing? Seems a little WIP, and perhaps only really used for internal stuff.

5. SequentialMonteCarlo.jl

6. SMC.jl (this one)

I’m sure there are more that I’m missing. I’m a bit of a SMC noob, so if anyone can fill me in on the current best package for doing what I’m trying to do, I would find it very helpful!

[baggepinnen]

Here’s something to get you started using LLPF, I hope it shouldn’t feel too involved. Where choices for parameters had to be made, I just picked some number out of a hat.

using LowLevelParticleFilters, Plots, Distributions, StaticArrays
# from julia epi-recipes


nx = 3   # Dimension of state
nu = 0   # Dimension of input
ny = 3   # Dimension of measurements
N = 500 # Number of particles
x0 = SA[990.0,10.0,0.0]
δt = 0.1
p = [0.05, 10.0, 0.25, δt]

const dg = MvNormal(ny,50.0)          # Measurement noise Distribution
const df = MvNormal(nx,1.0)          # Dynamics noise Distribution
const d0 = MvNormal(x0,0)   # Initial state Distribution

function rate_to_proportion(r, t)
    1-exp(-r*t)
end

function dynamics(x,u,p,t, noise=false)
    (S,I,R) = x
    (β,c,γ,δt) = p
    N = S+I+R
    ifrac = rate_to_proportion(β*c*I/N,δt)
    rfrac = rate_to_proportion(γ,δt)
    if noise
        infection=rand(Binomial(S,ifrac))
        recovery=rand(Binomial(I,rfrac))
    else
        infection = S*ifrac
        recovery = I*rfrac
    end
    SA[
        S-infection
        I+infection-recovery
        R+recovery
    ]
end;

function measurement_likelihood(x, u, y, p, t)
    logpdf(dg, x-y) # A simple linear measurement model with normal additive noise
end

measurement(x, u, p, t, noise=false) = max.(0, x + noise*rand(dg))

pf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0; p)

du = MvNormal(nu,1.0)         # Random input distribution for simulation
xs,u,y = simulate(pf,500,du; dynamics_noise=false, measurement_noise=true)  # Simulate data 500 time steps

sol = forward_trajectory(pf, u, y, p)

# plotly() # in case you want to zoom into the plot
plot(sol, xreal=xs)

[JasonPekos]

Thanks so much, this is really wonderful .