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| Xt - 1) p ( xt| X0 : t - 1, y1 : t)
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
p ( xt| y1 : t)tt
Considérez le modèle simple:
Xt= Xt - 1+ ηt,ηt∼ N( 0 , 1 )
X0∼ N( 0 , 1 )
Ouit= Xt+ εt,εt∼ N( 0 , 1 )
OuiXp ( Xt| Oui1, . . . , Yt)
Xt| Xt - 1∼ N( Xt - 1, 1 )
X0∼ N( 0 , 1)
Ouit| Xt∼ N(Xt, 1 )
Application de l'algorithme:
NX( je)0∼N( 0 , 1 )
X( i )1| X( i )0∼N( X( i )0, 1 )N
w~( i )t= ϕ ( yt; X( i )t, 1 )ϕ ( x ; μ , σ2)μσ2yt
wtXX( i )0 : t
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:
Un tutoriel utile est celui de Doucet et Johansen, voir ici .