# Best library for particle filtering from a generative model?

**URL:** <https://discourse.julialang.org/t/best-library-for-particle-filtering-from-a-generative-model/115070>\
**Category:** General Usage\
**Tags:** statistics, sciml, markov-decision\
**Created:** [June 1, 2024, 8:55pm UTC](https://discourse.julialang.org/t/best-library-for-particle-filtering-from-a-generative-model/115070 "2024-06-01T20:55:16Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![JasonPekos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jasonpekos/32/28438_2.png) [@JasonPekos](https://discourse.julialang.org/u/JasonPekos)\
**Post date:** [June 1, 2024, 8:55pm UTC](https://discourse.julialang.org/t/best-library-for-particle-filtering-from-a-generative-model/115070/1 "2024-06-01T20:55:16Z")

</div>

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:

> **[pomp, an R package for statistical inference using partially observed Markov...](https://kingaa.github.io/pomp/)**

or [particles](https://particles-sequential-monte-carlo-in-python.readthedocs.io/en/latest/notebooks/Bayes_estimation_ssm.html) 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:

```julia-auto
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:

```julia
# 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)](https://github.com/tlienart/SMC.jl?tab=readme-ov-file)

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!

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 2, 2024, 9:22am UTC](https://discourse.julialang.org/t/best-library-for-particle-filtering-from-a-generative-model/115070/2 "2024-06-02T09:22:18Z")

</div>

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.

```julia
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/5/553d02a26c0de2b5df3db4d099d6f0dbebaa2ce3.png)

---

<div class="post-metadata">

**Author:** ![JasonPekos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jasonpekos/32/28438_2.png) [@JasonPekos](https://discourse.julialang.org/u/JasonPekos)\
**Post date:** [June 2, 2024, 6:06pm UTC](https://discourse.julialang.org/t/best-library-for-particle-filtering-from-a-generative-model/115070/3 "2024-06-02T18:06:03Z")

</div>

Thanks so much, this is really wonderful 😃.