#'@export bootstrap_particle_filter
#'@title bootstrap_particle_filter
#'@description This function implements the bootstrap particle filter. It returns log-likelihood estimates, particles, and normalized weights.
#'@export
bootstrap_particle_filter <- function(observations, model, theta, algorithmic_parameters){
Nx <- algorithmic_parameters$Nx
resampling <- algorithmic_parameters$resampling
nobservations <- ncol(observations)
log_p_y_hat <- 0 #initialize estimate of p_theta(y_{1:nobservations}))
X_history = list() # successive sets of Nx particles column-wise
weight_history = list() # successive sets of normalized weights
incremental_ll = vector("numeric",nobservations) # successive incremental log-likelihood
X = matrix(NA,nrow = model$dimX, ncol = Nx) # matrix of Nx particles column-wise (most recent)
normW = rep(1/Nx, Nx) #vector of normalized weights
X <- model$rinitial(theta,Nx) #initial step 1
logW <- model$dobs(observations[,1],X,1,theta)
maxlogW <- max(logW)
W <- exp(logW - maxlogW)
normW <- W / sum(W)
incremental_ll[1] = log(mean(W)) + maxlogW
X_history[[1]] = X
weight_history[[1]] = normW
#iterate for n = 2, ... T
if (nobservations > 1){
for (t in 2:nobservations) {
ancestors <- resampling(normW) #sample the ancestors' indexes
X <- X[,ancestors]
if (is.null(dim(X))) X <- matrix(X, nrow = model$dimX)
X <- model$rtransition(X, t, theta)
logW <- model$dobs(observations[,t], X, t, theta)
maxlogW <- max(logW)
W <- exp(logW - maxlogW)
normW <- W / sum(W)
incremental_ll[t] = log(mean(W)) + maxlogW
X_history[[t]] = X
weight_history[[t]] = normW
}
}
return(list(log_p_y_hat = sum(incremental_ll), incremental_ll = incremental_ll,
X_history = X_history, weight_history = weight_history))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.