Basic Filtering Structure

[charlesknipp]

In the added example I propose a structure for generalized filtering, and additionally employ this method within PMMH given a random walk kernel.

Each filtering algorithm requires the user to construct 3 functions:

1. prior initialise for initializing the filtered states (whether it be particles or a Gaussian distribution)
2. predict which resamples the states and performs a one step ahead sampling step
3. update to evaluate the importance weight of the sample and return a marginal log-likelihood

Todo:

• fix adaptive resampling
• create stable resampling for Float32s
• improve type consistency for SSMs

Any feedback is appreciated, and considerable changes are more than welcome. This is really just a proof of concept, so feel free to propose any number of alterations.

[github-actions]

Remaining comments which cannot be posted as a review comment to avoid GitHub Rate Limit

JuliaFormatter

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Lines 231 to 234 in a5a2e05

	log_marginal = logpdf(
	Gaussian(zero(residual), Symmetric(S)),
	residual
	)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Line 249 in a5a2e05

resample_threshold(filter::BootstrapFilter) = filter.threshold*filter.N

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Lines 251 to 256 in a5a2e05

	function prior(
	rng::AbstractRNG,
	model::StateSpaceModel,
	filter::BootstrapFilter,
	extra
	)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Lines 258 to 261 in a5a2e05

	initial_states = map(
	x -> rand(rng, init_dist),
	1:filter.N
	)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Lines 267 to 274 in a5a2e05

	rng::AbstractRNG,
	particles::ParticleContainer,
	model::StateSpaceModel,
	filter::BootstrapFilter,
	step::Integer,
	extra;
	debug::Bool = false
	)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Line 286 in a5a2e05

idx = collect(1:filter.N)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Lines 290 to 291 in a5a2e05

	x -> SSMProblems.simulate(rng, model.dyn, step, x, extra),
	particles[idx]

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Lines 298 to 304 in a5a2e05

	particles::ParticleContainer,
	model::StateSpaceModel,
	observation,
	filter::BootstrapFilter,
	step::Integer,
	extra
	)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Line 307 in a5a2e05

collect(particles)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Lines 312 to 313 in a5a2e05

	ParticleContainer(particles.vals, particles.log_weights+log_marginals),
	logsumexp(log_marginals) - logsumexp(particles.log_weights)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

SSMProblems.jl/examples/particle-mcmc/simple-filters.jl

Line 315 in a5a2e05

end

[charlesknipp]

On Type Consistency

Currently, the interface contains type parameters for the state and observation objects; while this is ultimately good for forward simulation, we no longer have easy access to the element types. This is used mainly for preallocating particle weights, log evidence, etc and helps avoid unnecessary type conversions. At its worst, this raises errors when facing impossible type conversions.

Taking this to an extreme, I also think the user should keep the model element types consistent throughout the dynamics and observation process. This would allow us to employ a type structure like StateSpaceModel{T, LatentDynamics{Vector{T}}, ObservationProcess{Vector{T}}}, where Base.eltype(::StateSpaceModel) would return T. I would like to reiterate that keeping T consistent is the only way to avoid unnecessary/redundant type promotion.

[THargreaves]

Awesome stuff. I'll take a look through and make some comments.

Code style is Blue. I use the formatter baked into the Julia VSCode extension set up to format my code on every save. Works like a charm.

[THargreaves]

A few random thoughts:

• Should resampling take place in its own function? E.g. a step consists of resample, predict, update. For the tracking-in-clutter applications, it seems likely we're going to have to change step to a predict, associate, update loop, so doesn't feel as bad to add extra steps to the particle filter
• On the question of particle storage and ancestry (for naive smoothing), this may be too abstract but I wonder whether it is best left up to the storage object to determine how to save and store this. If the particle filter provides the new states, log weight increments and parent indexes, it can then be down to the storage object to store that as a linked list.
• As mentioned in my message about RTS, I wonder if rather than having a callback function, a callback struct with two methods, post_predict_store and post_update_store would be better. This would allow you to constantly update state in place rather than separating out into proposed_states and filtered_states.

I like the look of the "distribution-focused" version of the Kalman Filter. It might be worth leaving a comment clarifying what that is about for Frederic/Hong since I think we spoke more about that by ourselves. I'm apprehensive about using the name particles. state is much clear to me, where your state can be e.g. a Gaussian, a categorical distribution (HHM), or a weighted collection of point masses (particle collection).

[charlesknipp]

Particle Ancestry

I implemented the particle storage algorithm presented in (Murray, 2015). I include a demonstration at the bottom of this file which informally benchmarks the performance of the algorithm.

@THargreaves since you brought my attention to this algorithm, your comments would be much appreciated. I was aiming for elegance with this commit, and is thus far from optimal in terms of speed. Feel free to scrap anything remotely wasteful.

[THargreaves]

This looks great and closely matches the initial draft that I made a few months back.

I've pushed a small change that keeps track of the free indices using a stack. On my computer using the parameters in your test script, this reduced the median time from 48ms to 5ms and I assume will have better performance for sequences with fewer dead branches. I haven't tested it in full generality yet though but intuitively this feels it will be worth it most of the time.

Other than that, I think the particle storage code looks great.

I'm not sure about the predict and filter methods dispatching on AncestryContainer though. It would be much more general to have these implementation details abstracted away. Say by a store!(::AbstractParticleContainer, proposed_states idx) or something similar.

[charlesknipp]

I've pushed a small change that keeps track of the free indices using a stack. On my computer using the parameters in your test script, this reduced the median time from 48ms to 5ms...

@THargreaves legendary commit. 10x speed-up in less than 10 lines of changes.

I'm not sure about the predict and filter methods dispatching on AncestryContainer though. It would be much more general to have these implementation details abstracted away. Say by a store!(::AbstractParticleContainer, proposed_states idx) or something similar.

Yeah, it was purely to demonstrate a proof of concept. Now that you fixed the main drawback of the ancestry, I went ahead and added some key-word arguments to let initialise know whether to create an AncestryContainer or a ParticleContianer. It uses a nasty conditional, but I couldn't think of anything else.

I also consolidated my demonstrations to script.jl and moved resampling to a separate file.

[FredericWantiez]

We probably don't need the parent indices, that's up to the storage implementation to decide how they store ancestry paths:

Maybe something like that is enough ?

""" Store new generation in the container
"""
store!(pc::AbstractParticleContainer, new_generation)

""" Get last generation
"""
load(pc::AbstractParticleContainer)

""" Get (log)-weights from storage
"""
weights(pc::AbstractParticleContainer)
logweights(pc::AbstractParticleContainer)

[THargreaves]

We probably don't need the parent indices, that's up to the storage implementation to decide how they store ancestry paths

I might be misunderstanding you but I think we do, e.g. for smoothing.

If we just pass the particle container multiple vectors of states, it has no idea what the genealogy is so you can't perform naive smoothing on it by back-tracing ancestry.

Maybe something like that is enough ?

I'm a bit unsure on tis too. It feels like the filter is depending a bit too much on the implementation of the storage , which should ideally independent.

My instinct is for the filter to maintain a minimal collection of variables for it to run. I think this generally would just be the current state (represented here as a combination of x values and log weights). It would update these independently of the storage.

Then at each time step, it passes what it currently has to the storage object which can do what it wants with it. The key idea is that the filter should still work even in the extreme case that the storage throws everything away.

[THargreaves]

Actually I might be a bit confused here.

@charlesknipp, what was the intention of ParticleContainer? Is it a type of particle storage (one that just only remembers the current state) or just a representation of the current state?

If the latter, I don't think it and AncestoryContainer should be subtypes of the same type as they do different things.

[charlesknipp]

I currently use ParticleContainer as a means of storage to preserve the weighted nature of the sample at step t. Although I wonder if we could move ancestry storage to a callback, which would be very elegant if possible

[THargreaves]

Ah okay, I'm following now.

I didn't make this clear, but from my point of view store is the callback.

But rather than defining it as a simple function, it is tied to a storage container.

So the particle filter has a state::ParticleState(xs, log_ws) (currently called ParticleContainer but just to make the difference really clear) which it updates either in-place or by replacing with a new ParticleState and then this is passed to store! after each step to do what it pleases.

And store! can dispatch on SparseAncestoryStorage <: AbstractParticleStorage <: AbstractStorage or something like that, which is the Lawrence Murray algorithm you implemented.

[charlesknipp]

Exactly, right now store! is just a means of updating an AbstractStorage (or AbstractParticleContainer in my code).

I really like the idea you present with separating particle storage and particle state. Although, that would imply the need to store the ancestry indices in the particle state (which would be necessary for sparse ancestry storage). I'm not 100% sure of the details yet, but I think I can make this look pretty elegant.

[FredericWantiez]

tree::ParticleTree{T}, states::Vector{T}, ancestors::AbstractVector{Int64}