Simple particle filter demo

Preliminary stuff

1 Introduction

In this talk, we explore the connection between the brain’s ability to predict and higher cognition (the role of prediction in ‘higher’ cognition). In particular, we take a closer look at the much-discussed connection between prediction and mental imagery (Moulton & Kosslyn, Koziol, Pezzulo, Andy Clark, Gambi & Pickering).

We implement a simple agent, consisting of a sensor, an internal state and an algorithm allowing the agent to perform sensory inference. The internal state can encode prior knowledge in the form of beliefs about himself and the environment (the agent and the environment together constitute the world).

The primary goal of the agent is to perform sensory inference about his rotational velocity \(\omega\) in real time.

In order to perform real time inference, the agent performs particle filtering (sequential Monte Carlo, or SMC). We implement the simplest form of SMC, the bootstrap particle filter, or SMC with multinomial resampling.

The goal of the agent is simply to track his rotational velocity during a simple event. Events consist of sinusoidal acceleration to the left or the right:

\[ \frac{d\omega}{dt} = A \sin(2\pi\frac{t}{T}) \]

This type of acceleration reflects a natural head movement to one side, and within a frequency band of \(0.1 - 100\) [check this], the SCC provides measurement signal that are proportional to the velocity:

\[ \omega = A\frac{T}{2\pi}\Big(1-\cos(2\pi \frac{t}{T})\Big) \]

We choose a frequency \(\mathcal{f}\) of 0.5 Hz. The agent’s sensor provides a sequence of noisy measurements \(y_{t}\) of his angular velocity \(\omega^{true}_{t}\):

\[ y_{t} \sim \mathcal{N}(\omega^{true}_{t}, \sigma^{2}_{sensor}) \]

where \(\omega^{true}_{t}\) is the actual angular velocity, and \(\sigma^{2}_{sensor}\) is the noise of the sensor (physical characteristic of the sensor).

Figure 1 shows the acceleration, velocity, resulting angular deflection, and the agent’s noisy sensory measurements:

set.seed(8529573L)
data <- generate_data(T = 2, amplitude = 20, sensor_sd = 1.7) %>% 
    t() %>% 
    data.frame()
plot_trajectories(data = data)

For all simulation, we use a sensor standard deviation \(\sigma^{2}_{sensor}\) of \(1.7^\circ\).

1.1 Particle filter implementation

We will go step-by-step through the basic SMC algorithm. The agents goal is to infer the current angular velocity \(\omega_{t}\), which is a latent variable, given noisy sensory measurements of \(\omega_{t}\). The agent thus requires a generative model of its measurements, which in this case is given by:

\[ y_{t} \sim \mathcal{N}(\omega_{t}, \sigma^2_y) \]

This is the agent’s measurement model:

• \(\sigma^2_y\) is the agent’s model of the sensor noise - this is distinct from the actual (true) sensor noise \(\sigma^2_{sensor}\).

In addition, the agent must have a model of the evolution of its state variable through time, i.e. a model of the dynamics of the latent state \(\omega\). This is given by:

\[ \omega_{t} \sim \mathcal{N}(\omega_{t-1} + A \sin(\frac{2\pi t}{T}), \sigma^2_{process}) \]

where \(A \sin(\frac{2\pi t}{T})\) represents a ‘known’ control input and \(\sigma^2_{process}\) is the process standard deviation. This term encodes the agent’s certainty in his process model. The control input is assumed to have been learned by the agent, and represent knowledge of the sensory consequences of active head movements. This represents the concept typically known as an ‘efference copy’, although in our model this is merely a statistical model, of the sensory effects of motor action, and is thus equivalent to a ‘known’ covariate that the agent can use to predict the latent state.

1.2 Particle filter step-by-step (during real-time interaction with the environment)

prior to filtering: The agent chooses parameters for the filtering algorithm:

sdx <- 0.5 # process noise
sdy <- 2 # sensor noise
A <- 2 # motion amplitude

The agent also decides which process model to use:

# process model (transition model)
f <- function(x, Time, A, sd, N) {
    rnorm(n = N, mean =  x[, ] + A * sin(2*pi/Time), sd = sd) 
}

In addition, the agent decides how many resources to allocate to this particular problem, in the form of the number of particles (samples) to use.

at \(t_0\): The agent starts by making an initial guess about \(\omega\). This means that the agent draws a sample from the initial distribution; in our model, we assume that the agent starts in an unbiased state with a high degree of certainty:

In the following code, the state variable \(\omega_t\) is denoted as \(x_t\).

# x is drawn from the prior distribution with 
# expected value 0 and a small standard deviation 0.1)
# The agent has decided to use 100 particles
N <- 100
x_init <- rnorm(N, 0, 0.1)

The agent thus represents his belief about the initial velocity as a collection of 100 samples from the prior distribution.

The agent then predicts the state at the next time increment \(t_1\) by applying the process model to each of the \(N\) particles. The process model reflects the agents knowledge about Newtonian dynamics:

Time <- 20
# N <- 100
x <- matrix(rep(0, N*Time), nrow = N, ncol = Time)  #matrix(nrow =  N, ncol = Time)
weights <- matrix(rep(0, N*Time), nrow = N, ncol = Time)
xmean <- rep(0, Time) # mean of estimated state
xvar <- rep(0, Time) # variance of estimated state
x[, 1] <- rnorm(n = N, mean =  x_init + A * sin(2*pi/Time), sdx)
# x[, 1] <- f(x_init, Time, A, sd = sdx, N)

In MCMC terminology, the predicted states represent samples from a proposal distribution.

at \(t_1\): There are sensory measurements available, so the agent can use these in order to update his internal estimate of his state. To do so, he evaluates the likelihood of each particle, i.e. he evaluated the probability of observing the measurement at time \(t_1\) given the predicted velocity:

weights[, 1] <- dnorm(data$observations[1], x[, 1], sdy)

The probabilities are then used as importance weights for each particle. Intuitively, this means that the agent selects those particles that predicted the data well, and discards those particles that predicted the data poorly. This is achieved by a weighted sampling with replacement from the set of particles (multinomial resampling). The weights are first normalized, so that they represent a probability distribution.

x[, 1] <- sample(x[, 1], replace = TRUE, size = N, prob = weights[, 1])
xmean[1] <- mean(x[, 1])
xvar <- var(x[, 1])

By applying this simple procedure recursively, the agent can track his latent velocity, given the sequence of noisy sensory measurements. At any time step, the agent thus performs the following sequence of computations:

1. prediction (propagation) step: predict the state 1 step ahead
2. update step: obtain importance weights from sensory measurements, and select particles for ‘survival’
3. compute summary statistics, e.g. mean and variance of estimated state

Figure 2 shows a visualization of the computations performed at each time step. After the prediction step, the set of particles represent the new belief before ‘seeing’ the new evidence. All particles thus have equal weights, given \(\frac{1}{N}\). Subsequently, the probability of observing the evidence is evaluated for each particle, giving new weight. Each particle’s size reflects its importance, i.e. how well it predicted the new measurement. The final step consists of sampling 100 new particles with replacement from the set of particles, with a probability of being chosen given by their importance weights. (This is the multinomial resampling part). Particles that did not predict the data well may therefore be eliminated.

The purple line shows the estimated mean of the distribution being ‘shifted’ during a single time increment, and the black rectangle show the observations; at the first step, before seeing the new data, the observation is at the previous value. At the updating step (step 2), a new measurement is available. The closer to the measurement, the larger the diameter of the particle.

Warning: Removed 1 rows containing missing values (geom_point).

Figure 2 demonstrates how, during real-time interaction with the environment, 1-step ahead prediction is performed continuously and is a vital ingredient in (normal) sensory inference.

1.3 Predicting k steps ahead (what happens when data are unavailable)

Next, we consider an interaction episode during which sensory data temporarily become unavailable.

In the visual system, this is fairly common: lights out, eyes closed, fog, etc. What are analogies in the the vestibular system?

Figure 3 shows the agent enganged in real-time sensory inference during a 2 second motion event. The dotted line shows the true velocity, the black rectangles are the sensory measurements, the solid line show the agents estimated head velocity, and the bands surrounding the solid line reflect the agents \(95\%\) credible interval.

Due to the process of predicting and updating, the confidence region remains approximately constant throughout the event. Predicting increases uncertainty, updating reduces it again.

Figure 4 shows the same sensory event, but after 1 seconds, there is sensor failure, and the agent no longer receives any measurements. Forced to rely on his internal model, the agent is only able to predict the next steps, without being able to update. This is implemented by letting the agent monitor whether data are available, and then performing weighting according to the likelihood if they are. IN the case of missing data, the agent simply assigns a uniform probability distribution to the weights. Resampling is subsequently performed, but if the data are missing, this has no effect, because the resampled set of particles are simply drawn with equal probabilities from the old set. The particles therefore do not move torward the target distribution, and the variance of the samples increases over time. Note that the estimate remains unbiased, because the agent is using a ‘good’ generative model, i.e. the agent is correctly modelling the cause of the true velocity. This is illustrated by the fact that the estimated line lies ‘on top of’ the true velocity. (In fact, if the agent is using a good process model, the recovered estimate is actually better if he is not forced to incorporate noisy measurements; maybe the figure should not show this?).

Many would claim that this is already imagery, of guided hallucination…

run_particle_filter <- function(data, N, Time, x_init, sdx_init,
                                params, resample = TRUE, rs_thresh = 0.5) {
    
    
    sdc <- params$sdc
    sdx <- params$sdx
    sdy <- params$sdy
    fun_x <- params$funx_x
    fun_c <- params$fun_c
    A <- params$A
    
    x <- matrix(rep(0, N*Time), nrow = N, ncol = Time)  #matrix(nrow =  N, ncol = Time)
    weights <- matrix(rep(0, N*Time), nrow = N, ncol = Time) #matrix(nrow =  N, ncol = Time)
    loglik <- rep(0, Time)
    
    x[, 1] <- rnorm(N, x_init, sdx_init)
    weights[, 1] <- dnorm(data$observations[1], x[, 1], sdy)
    weights[, 1] <- weights[, 1]/sum(weights[, 1])
    x[, 1] <- sample(x[, 1], replace = TRUE, size = N, prob = weights[, 1])
    
    
    for (t in seq(2, Time)) {
        x[, t] <- f(x, t, Time, A, sd = sdx, N)
        
        if (!is.na(data$observations[t])) {
            weights[, t] <- dnorm(data$observations[t], x[, t], sdy)
        } else {
            weights[, t] <- 1/N
        }
        
        if (resample) {
                if (neff(weights[, t]) < rs_thresh * N) {
                    x[, t] <- sample(x[, t], replace = TRUE, size = N, prob = weights[, t]) 
            }
        }
        
        # TODO: do we need to reset weights after resampling?
        # weights[, t] <- 1/N
    }
    
    out <- list(x = x, weights = weights, loglik = loglik,
                x_means = apply(x, 2, mean),
                x_medians = apply(x, 2, median),
                x_quantiles = apply(x, 2, 
                                    function(x) quantile(x, 
                                                         probs = c(0.025, 
                                                                   0.975))),
                Time = Time,
                N = N)
}

f <- function(x, t, Time, A, sd, N) {
    rnorm(n = N, mean =  x[, t-1] + A * sin(2*pi*(t-1)/Time), sd = sd) 
}
f1 <- function(x, t, Time, A, sd, N) {
    rnorm(n = N, mean =  x[, t-1], sd = sd) 
}
f2 <- function(x, t, Time, A, sd, N) {
    # A * sin(2*pi*(t-1)/Time)
    rnorm(n = N, mean =  A * sin(2*pi*(t-1)/Time), sd = sd) 
}

params <- list(sdx = 1.2, sdy = 4.0, A = 2.0, fun_x = f)
out <- run_particle_filter(data = data, N = 100, Time = length(data$observations),
                           x_init = 0, sdx_init = 0.5,
                           params, resample = TRUE, rs_thresh = 0.5)
plot_filtering_estimates(out, data = data)

params <- list(sdx = 1.2, sdy = 4.0, A = 2.0, fun_x = f)
data_missing <- data
data_missing$observations[10:length(data_missing$observations)] <- NA
out <- run_particle_filter(data = data_missing, N = 100, Time = length(data$observations),
                           x_init = 0, sdx_init = 0.5,
                           params, resample = TRUE, rs_thresh = 0.5)
plot_filtering_estimates(out, data = data_missing, predict = TRUE)

1.4 Running the particle filter with a misspecified process model

The next two Figures show what can happen when running the particle filter using a ‘wrong’ process model. Here, the agent believes that the direction of motion is in the opposite direction to the actual motion, leading to a misinterpretation.

params <- list(sdx = 1.2, sdy = 4.0, A = -2, fun_x = f)
out <- run_particle_filter(data = data, N = 100, Time = length(data$observations),
                           x_init = 0, sdx_init = 0.5,
                           params, resample = TRUE, rs_thresh = 0.5)
plot_filtering_estimates(out, data = data)

In the first Figure, the agent has a high degree of certainy in his process model: \(\sigma_{process}\) = 2.8.

params <- list(sdx = 2.2, sdy = 4.0, A = -2, fun_x = f)
out <- run_particle_filter(data = data, N = 100, Time = length(data$observations),
                           x_init = 0, sdx_init = 0.5,
                           params, resample = TRUE, rs_thresh = 0.5)
plot_filtering_estimates(out, data = data)

In the second Figure, the agent has a lower degree of certainy in his process model: \(\sigma_{process}\) = 2.8.

Thus, using the same process model, with the same sensory input, the agent may interpret the data differently, depending on how certain he is in his own predictions (due to his process model).

Question: But is this already enough for imagery/SPT?

I think not, because in all the above scenarios, the agent was still engaged in real-time interaction with the environment, i.e. the agent’s goal was to explain the sensory. IN all cases, he wanted to uncover the causes of his sensory input, and was able to allow for missing data, or used a wrong process model.

1.5 Running the particle filter offline (without resampling)

Now comes the important next step: in order to perform SPT, etc., the agent must actively disengage from his real-time interaction with the environment, otherwise he would perform compensatory eye movements, etc. But what can the agent do in order to go offline? Of course, assuming the agent resuses his sensory circuit for ‘cognition’.

In addition to the previous abilities:

• ability to predict (1 step ahead)
• ability to deal with missing data (predict multiple (k) steps ahead)
• ability to apply a process model (control signal)

the agent must perform an active intervention in his probabilistic model; our proposal is that he can change his resampling strategy - he can choose not to resample. By not resampling he can ignore the incoming data.

An advantage of this strategy is that he is running a model, and the sensory data are still coming in. Therefore, he can still monitor the likelihood of his model, i.e. he can do a reality check (hypothesis test) to make sure that his ‘fantasy’ model does not explain reality, i.e. it is not real.

The next Figure shows the agent’s sensory input whilst remaining stationary.

set.seed(8529573L)
data_offline <- generate_data(T = 2, amplitude = 0, sensor_sd = 1.7) %>% 
    t() %>% 
    data.frame()
plot_trajectories(data = data_offline)

The last two Figures show the agent shows the agent engaged in ‘imagining’ a rightward motion, with amplitude \(A\) = -2, using two different degrees of certainty. In the first, the agents

params <- list(sdx = 0.8, sdy = 4.0, A = 5, fun_x = f)
out <- run_particle_filter(data = data_offline, N = 100, Time = length(data$observations),
                           x_init = 0, sdx_init = 0.5,
                           params, resample = FALSE, rs_thresh = 0.5)
plot_filtering_estimates(out, data = data_offline)

In the first, \(\sigma_{process}\) = 2.8; therefore, the agent is very certain of his belief, i.e. he is experiencing vivid imagery.

params <- list(sdx = 2.8, sdy = 4.0, A = -2, fun_x = f)
out <- run_particle_filter(data = data_offline, N = 100, Time = length(data$observations),
                           x_init = 0, sdx_init = 0.5,
                           params, resample = FALSE, rs_thresh = 0.5)
plot_filtering_estimates(out, data = data_offline)

In the second, \(\sigma_{process}\) = 2.8; therefore, the agent is very uncertain of his belief, i.e. his imagery is not very vivid.

Note, however, that in both cases, the agent is not interested in ‘explaining’ the sensory data, i.e. his imagined motion is not constrained by the measurements.