# Particle Filtering in Gen
## Tutorial: Particle Filtering in Gen
Sequential Monte Carlo (SMC) methods such as particle filtering iteratively solve a *sequence of inference problems* using techniques based on importance sampling and in some cases MCMC. The solution to each problem in the sequence is represented as a collection of samples or *particles*. The particles for each problem are based on extending or adjusting the particles for the previous problem in the sequence.
The sequence of inference problems that are solved often arise naturally from observations that arrive incrementally, as in particle filtering. The problems can also be constructed instrumentally to facilitate inference, as in annealed importance sampling \[3]. This tutorial shows how to use Gen to implement a particle filter for a tracking problem that uses “rejuvenation” MCMC moves. Specifically, we will address the “bearings only tracking” problem described in \[4].
\[1] Doucet, Arnaud, Nando De Freitas, and Neil Gordon. “An introduction to sequential Monte Carlo methods.” Sequential Monte Carlo methods in practice. Springer, New York, NY, 2001. 3-14.
\[2] Del Moral, Pierre, Arnaud Doucet, and Ajay Jasra. “Sequential monte carlo samplers.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68.3 (2006): 411-436.
\[3] Neal, Radford M. “Annealed importance sampling.” Statistics and computing 11.2 (2001): 125-139.
\[4] Gilks, Walter R., and Carlo Berzuini. “Following a moving target—Monte Carlo inference for dynamic Bayesian models.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63.1 (2001): 127-146. [PDF](http://www.mathcs.emory.edu/~whalen/Papers/BNs/MonteCarlo-DBNs.pdf)
### Outline
**Section 1**: [Implementing the generative model](https://probcomp.github.io/Gen/particle-filtering-in-gen/Particle%20Filtering%20in%20Gen#basic-model)
**Section 2**: [Implementing a basic particle filter](https://probcomp.github.io/Gen/particle-filtering-in-gen/Particle%20Filtering%20in%20Gen#basic-pf)
**Section 3**: [Adding rejuvenation moves](https://probcomp.github.io/Gen/particle-filtering-in-gen/Particle%20Filtering%20in%20Gen#rejuv)
**Section 4**: [Using the unfold combinator to improve performance](https://probcomp.github.io/Gen/particle-filtering-in-gen/Particle%20Filtering%20in%20Gen#unfold)
```
using Gen
using PyPlot
```
### 1. Implementing the generative model
We will implement a generative model for the movement of a point in the x-y plane and bearing measurements of the location of this point relative to the origin over time.
We assume that we know the approximate initial position and velocity of the point. We assume the point’s x- and y- velocity are subject to random perturbations drawn from some normal distribution with a known variance. Each bearing measurement consists of the angle of the point being tracked relative to the positive x-axis.
We write the generative model as a generative function below. The function first samples the initial state of the point from a prior distribution, and then generates `T` successive states in a `for` loop. The argument to the model is the number of states not including the initial state.
```
bearing(x, y) = atan(y, x)
@gen function model(T::Int)
measurement_noise = 0.005
velocity_var = (1.0/1e6)
xs = Vector{Float64}(undef, T+1)
ys = Vector{Float64}(undef, T+1)
# prior on initial x-coordinate
x = @trace(normal(0.01, 0.01), :x0)
# prior on initial y-coordinate
y = @trace(normal(0.95, 0.01), :y0)
# prior on x-component of initial velocity
vx = @trace(normal(0.002, 0.01), :vx0)
# prior on y-component of initial velocity
vy = @trace(normal(-0.013, 0.01), :vy0)
# initial bearing measurement
@trace(normal(bearing(x, y), measurement_noise), :z0)
# record position
xs[1] = x
ys[1] = y
# generate successive states and measurements
for t=1:T
# update the state of the point
vx = @trace(normal(vx, sqrt(velocity_var)), (:vx, t))
vy = @trace(normal(vy, sqrt(velocity_var)), (:vy, t))
x += vx
y += vy
# bearing measurement
@trace(normal(bearing(x, y), measurement_noise), (:z, t))
# record position
xs[t+1] = x
ys[t+1] = y
end
# return the sequence of positions
return (xs, ys)
end;
```
We generate a data set of positions, and observed bearings, by sampling from this model, with `T=50`:
```
import Random
Random.seed!(4)
# generate trace with specific initial conditions
T = 50
constraints = Gen.choicemap((:x0, 0.01), (:y0, 0.95), (:vx0, 0.002), (:vy0, -0.013))
(trace, _) = Gen.generate(model, (T,), constraints)
# extract the observed data (zs) from the trace
choices = Gen.get_choices(trace)
zs = Vector{Float64}(undef, T+1)
zs[1] = choices[:z0]
for t=1:T
zs[t+1] = choices[(:z, t)]
end
```
We next write a visualization for traces of this model below. It shows the positions and dots and the observed bearings as lines from the origin:
```
function render(trace; show_data=true, max_T=get_args(trace)[1])
(T,) = Gen.get_args(trace)
choices = Gen.get_choices(trace)
(xs, ys) = Gen.get_retval(trace)
zs = Vector{Float64}(undef, T+1)
zs[1] = choices[:z0]
for t=1:T
zs[t+1] = choices[(:z, t)]
end
scatter(xs[1:max_T+1], ys[1:max_T+1], s=5)
if show_data
for z in zs[1:max_T+1]
dx = cos(z) * 0.5
dy = sin(z) * 0.5
plot([0., dx], [0., dy], color="red", alpha=0.3)
end
end
end;
```
We visualize the synthetic trace below:
```
render(trace)
```
### 2. Implementing a basic particle filter
In Gen, a **particle is represented as a trace** and the particle filter state contains a weighted collection of traces. Below we define an inference program that runs a particle filter on an observed data set of bearings (`zs`). We use `num_particles` particles internally, and then we return a sample of `num_samples` traces from the weighted collection that the particle filter produces.
Gen provides methods for initializing and updating the state of a particle filter, documented in [Particle Filtering](https://probcomp.github.io/Gen/dev/ref/inference/#Particle-Filtering-1).
* `Gen.initialize_particle_filter`
* `Gen.particle_filter_step!`
Both of these methods can used either with the default proposal or a custom proposal. In this tutorial, we will use the default proposal. There is also a method that resamples particles based on their weights, which serves to redistribute the particles to more promising parts of the latent space.
* `Gen.maybe_resample!`
Gen also provides a method for sampling a collection of unweighted traces from the current weighted collection in the particle filter state:
* `Gen.sample_unweighted_traces`
```
function particle_filter(num_particles::Int, zs::Vector{Float64}, num_samples::Int)
# construct initial observations
init_obs = Gen.choicemap((:z0, zs[1]))
state = Gen.initialize_particle_filter(model, (0,), init_obs, num_particles)
# steps
for t=1:length(zs)-1
Gen.maybe_resample!(state, ess_threshold=num_particles/2)
obs = Gen.choicemap(((:z, t), zs[t+1]))
Gen.particle_filter_step!(state, (t,), (UnknownChange(),), obs)
end
# return a sample of unweighted traces from the weighted collection
return Gen.sample_unweighted_traces(state, num_samples)
end;
```
The initial state is obtained by providing the following to `initialize_particle_filter`:
* The generative function for the generative model (`model`)
* The initial arguments to the generative function.
* The initial observations, expressed as a map from choice address to values (`init_obs`).
* The number of particles.
At each step, we resample from the collection of traces (`maybe_resample!`) and then we introduce one additional bearing measurement by calling `particle_filter_step!` on the state. We pass the following arguments to `particle_filter_step!`:
* The state (it will be mutated)
* The new arguments to the generative function for this step. In our case, this is the number of measurements beyond the first measurement.
* The [argdiff](https://probcomp.github.io/Gen/dev/ref/gfi/#Argdiffs-1) value, which provides detailed information about the change to the arguments between the previous step and this step. We will revisit this value later. For now, we indicat ethat we do not know how the `T::Int` argument will change with each step.
* The new observations associated with the new step. In our case, this just contains the latest measurement.
We run this particle filter with 5000 particles, and return a sample of 100 particles. This will take 30-60 seconds. We will see one way of speeding up the particle filter in a later section.
```
@time pf_traces = particle_filter(5000, zs, 200);
```
```
36.481359 seconds (134.90 M allocations: 5.846 GiB, 47.35% gc time)
```
To render these traces, we first define a function that overlays many renderings:
```
function overlay(renderer, traces; same_data=true, args...)
renderer(traces[1], show_data=true, args...)
for i=2:length(traces)
renderer(traces[i], show_data=!same_data, args...)
end
end;
```
We then render the traces from the particle filter:
```
overlay(render, pf_traces, same_data=true)
```
Notice that as during the period of denser bearing measurements, the trajectories tend to turn so that the heading is more parallel to the bearing vector. An alternative explanation is that the point maintained a constant heading, but just slowed down significantly. It is interesting to see that the inferences favor the “turning explanation” over the “slowing down explanation”.
---
# Agent Instructions: Querying This Documentation
If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.
Perform an HTTP GET request on the current page URL with the `ask` query parameter:
```
GET https://alfredo-reyes-montero.gitbook.io/gen/particle-filtering-in-gen.md?ask=
```
The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.
Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.