Filtre bootstrap / algorithme de filtre à particules (compréhension)

12

J'ai vraiment un manque de compréhension du fonctionnement du filtre bootstrap. Je connais à peu près les concepts mais je n'arrive pas à saisir certains détails. Cette question est pour moi d'éliminer l'encombrement. Ici, je vais utiliser cet algorithme de filtre populaire à partir d'une référence par doucet (jusqu'à présent, je pense que c'est la référence la plus simple). Permettez-moi d'abord de vous dire que mon problème est de comprendre quelles distributions sont connues et lesquelles sont inconnues.

filtre bootstrap

Ce sont mes questions:

  1. En 2), quelle est la distribution ? Cette distribution est-elle connue ? Connaissons-nous cette distribution pour tous les ? Si oui, mais que faire si nous ne pouvons pas en échantillonner? C'est drôle qu'ils appellent cette étape d'échantillonnage d'importance, mais je ne vois aucune distribution de proposition.tp(xt|xt1(i))t
  2. Également en 2) est une distribution connue ? "Normaliser les poids d'importance signifie ? Que signifie le tilde sur et ? Cela signifie-t-il quelque chose comme non rééchantillonné ou non normalisé respectivement?w ( i ) t = ˜ w ( i ) tp(yt|x~t(i)) xwwt(i)=w~t(i)i=1Nw~t(i)xw
  3. J'apprécierais que quelqu'un puisse donner un exemple de jouet simple en utilisant des distributions bien connues pour utiliser ce filtre d'amorçage. L'objectif final du filtre d'amorçage n'est pas clair pour moi.
tintinthong
la source

Réponses:

10
  1. Il s'agit de la densité de transition de l'état ( ), qui fait partie de votre modèle et est donc connue. Vous devez en échantillonner dans l'algorithme de base, mais des approximations sont possibles. p ( x t | x t - 1 ) est la distribution de la proposition dans ce cas. Il est utilisé car la distribution p ( x t | x 0 : t - 1 , y 1 : t ) n'est généralement pas traitable.xtp(xt|xt1) p(xt|x0:t1,y1:t)

  2. Oui, c'est la densité d'observation, qui fait également partie du modèle, et donc connue. Oui, c'est ce que signifie la normalisation. Le tilde est utilisé pour signifier quelque chose comme "préliminaire": est x avant le rééchantillonnage, et ˜ w est w avant la renormalisation. Je suppose que cela se fait de cette façon afin que la notation corresponde aux variantes de l'algorithme qui n'ont pas d'étape de rééchantillonnage (c'est-à-dire que x est toujours l'estimation finale).x~xw~wx

  3. p(xt|y1:t)tt

Considérez le modèle simple:

Xt=Xt1+ηt,ηtN(0,1)
X0N(0,1)
Yt=Xt+εt,εtN(0,1)

YXp(Xt|Y1,...,Yt)

Xt|Xt1N(Xt1,1)
X0N(0,1)
Yt|XtN(Xt,1)

Application de l'algorithme:

  1. NX0(i)N(0,1)

  2. X1(i)|X0(i)N(X0(i),1)N

    w~t(i)=ϕ(yt;xt(i),1)ϕ(x;μ,σ2)μσ2yt

  3. wtxx0:t(i)

Revenez à l'étape 2, en avançant avec la version rééchantillonnée des particules, jusqu'à ce que nous ayons traité toute la série.

Une implémentation en R suit:

# 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)

Le graphique résultant:entrez la description de l'image ici

Un tutoriel utile est celui de Doucet et Johansen, voir ici .

Chris Haug
la source
X1(i)|X0(i)N(0,1)X1(i)|X0(i)N(X0(i),1)
tintinthong
C'est exact, j'ai corrigé la faute de frappe
Chris Haug
Les chemins n'ont pas besoin d'être rééchantillonnés, n'est-ce pas ?? À partir d'autres publications, il n'est pas nécessaire d'échantillonner les chemins. J'ai juste besoin d'échantillonner les particules à chaque pas de temps. Je me demandais s'il y avait une raison de rééchantillonner les chemins
tintinthong