Algoritmo de filtro de arranque / filtro de partículas (comprensión)

Realmente no entiendo cómo funciona el filtro de arranque. Conozco los conceptos más o menos, pero no entiendo ciertos detalles. Esta pregunta es para mí para aclarar el desorden. Aquí usaré este algoritmo de filtro popular de una referencia de doucet (hasta ahora creo que esta es la referencia más fácil). Permítanme decirles primero que mi problema es comprender qué distribuciones se conocen y cuáles se desconocen.

Estas son mis preguntas:

1. En 2), ¿cuál es la distribución ? ¿Se conoce esta distribución ? ¿Conocemos esta distribución para todas las ? Si es así, ¿pero qué pasa si no podemos probarlo? Es curioso que llamen a este paso de muestreo importante, pero no veo distribución de propuestas.$p(x_t|x^{(i)}_{t-1})$$t$
2. También en 2) ¿es una distribución conocida ? "Normalizar pesos de importancia significa ? ¿Qué significa tilde en y ? ¿Significa algo como no remuestreado o no normalizado respectivamente?w ( i ) t = ˜ w ( i ) $p(y_t|\tilde{x}^{(i)}_{t})$ x$w^{(i)}_{t}=\frac{\tilde{w}^{(i)}_{t}}{\sum_{i=1}^{N}\tilde{w}^{(i)}_{t}}$$x$$w$
3. Agradecería si alguien pudiera dar un ejemplo simple de juguete usando distribuciones bien conocidas para usar este filtro de arranque. El objetivo final del filtro de arranque no está claro para mí.

1. Esa es la densidad de transición del estado ( ), que es parte de su modelo y, por lo tanto, conocida. Es necesario tomar muestras de él en el algoritmo básico, pero son posibles aproximaciones. es la distribución de la propuesta en este caso. Se utiliza porque la distribución generalmente no es manejable.$x_t$$p(x_t|x_{t-1})$ $p(x_t|x_{0:t-1},y_{1:t})$

2. Sí, esa es la densidad de observación, que también es parte del modelo y, por lo tanto, conocida. Sí, eso es lo que significa la normalización. La tilde se usa para significar algo así como "preliminar": es antes del remuestreo, y es antes de la renormalización. Supongo que se hace de esta manera para que la notación coincida entre las variantes del algoritmo que no tienen un paso de remuestreo (es decir, es siempre la estimación final).$\tilde{x}$$x$$\tilde{w}$$w$$x$

3. El objetivo final del filtro de arranque es estimar la secuencia de distribuciones condicionales (el estado no observable en , dadas todas las observaciones hasta ).$p(x_t|y_{1:t})$$t$$t$

Considere el modelo simple:

X 0 ∼ N ( 0 , 1 ) Y t = X t + ε t

Esta es una caminata aleatoria observada con ruido (solo observas , no ). Puede calcular exactamente con el filtro de Kalman, pero usaremos el filtro de arranque a petición suya. Podemos reformular el modelo en términos de la distribución de transición de estado, la distribución de estado inicial y la distribución de observación (en ese orden), que es más útil para el filtro de partículas:$Y$$X$$p(X_t|Y_1, ..., Y_t)$

Aplicando el algoritmo:

1. Inicializacion. Generamos partículas (independientemente) de acuerdo con .$N$$X_0^{(i)} \sim N(0,1)$

2. Simulamos cada partícula hacia adelante independientemente generando , para cada .$X_1^{(i)} | X_0^{(i)} \sim N(X_0^{(i)},1)$$N$

   Luego calculamos la probabilidad , donde es el densidad normal con media y varianza (nuestra densidad de observación). Queremos dar más peso a las partículas que tienen más probabilidades de producir la observación que registramos. Normalizamos estos pesos para que sumen 1.$\tilde{w}_t^{(i)} = \phi(y_t; x_t^{(i)},1)$$\phi(x; \mu, \sigma^2)$$\mu$$\sigma^2$$y_t$

3. Volvemos a muestrear las partículas de acuerdo con estos pesos . Tenga en cuenta que una partícula es una ruta completa de (es decir, no solo muestrea el último punto, es todo, lo que denotan como ).$w_t$$x$$x_{0:t}^{(i)}$

Regrese al paso 2, avanzando con la versión muestreada de las partículas, hasta que hayamos procesado toda la serie.

Sigue una implementación en R:

# Simulate some fake data
set.seed(123)

tau <- 100
x <- cumsum(rnorm(tau))
y <- x + rnorm(tau)

# Begin particle filter
N <- 1000
x.pf <- matrix(rep(NA,(tau+1)*N),nrow=tau+1)

# 1. Initialize
x.pf[1, ] <- rnorm(N)
m <- rep(NA,tau)
for (t in 2:(tau+1)) {
  # 2. Importance sampling step
  x.pf[t, ] <- x.pf[t-1,] + rnorm(N)

  #Likelihood
  w.tilde <- dnorm(y[t-1], mean=x.pf[t, ])

  #Normalize
  w <- w.tilde/sum(w.tilde)

  # NOTE: This step isn't part of your description of the algorithm, but I'm going to compute the mean
  # of the particle distribution here to compare with the Kalman filter later. Note that this is done BEFORE resampling
  m[t-1] <- sum(w*x.pf[t,])

  # 3. Resampling step
  s <- sample(1:N, size=N, replace=TRUE, prob=w)

  # Note: resample WHOLE path, not just x.pf[t, ]
  x.pf <- x.pf[, s]
}

plot(x)
lines(m,col="red")

# Let's do the Kalman filter to compare
library(dlm)
lines(dropFirst(dlmFilter(y, dlmModPoly(order=1))$m), col="blue")

legend("topleft", legend = c("Actual x", "Particle filter (mean)", "Kalman filter"), col=c("black","red","blue"), lwd=1)

El gráfico resultante:

Un tutorial útil es el de Doucet y Johansen, ver aquí .