R/particle_filter_storeall.R

In pierrejacob/CoupledPF: Coupling of Particle Filters

#'@rdname particle_filter_storeall
#'@title Particle Filter
#'@description Runs a particle filter, storing all the generated variables
#'@param nparticles number of particles
#'@param model a list representing a model, for instance as given by \code{\link{get_ar}}.
#'@param theta a parameter to give to the model functions
#'@param observations a matrix of observations of size datalength x dimension(observation)
#'@param randomness an object storing all the variables useful to initialize and propagate the particles,
#' as generated by a model$generate_randomness function 
#'@return A list with entries ll (loglikelihood estimator), xhistory, ahistory, whistory
#' (lists containing all the generated x-locations, ancestors and weights)
#'@export
particle_filter_storeall <- function(nparticles, model, theta, observations, randomness){
  datalength <- nrow(observations)
  precomputed <- try(model$precompute(theta))
  if (inherits(precomputed, "try-error")){
    precomputed <- NULL
  }
  # initialization
  xparticles <- model$rinit(nparticles, theta, randomness)
  normweights <- rep(1/nparticles, nparticles)
  ll <- 0
  #
  xhistory <- rep(list(matrix(ncol = nparticles, nrow = nrow(xparticles))), datalength + 1)
  whistory <- rep(list(rep(0, nparticles)), datalength + 1)
  ahistory <- rep(list(rep(0, nparticles)), datalength)
  xhistory[[1]] <- xparticles
  whistory[[1]] <- normweights
  # step t > 1
  for (time in 1:datalength){
    ancestors <- systematic_resampling_n(normweights, nparticles, runif(1))
    ahistory[[time]] <- ancestors
    xparticles <- xparticles[,ancestors]
    if (is.null(dim(xparticles))) xparticles <- matrix(xparticles, nrow = model$dimension)
    xparticles <- model$rtransition(xparticles, theta, time, randomness, precomputed)
    logw <- model$dmeasurement(xparticles, theta, observations[time,], precomputed)
    maxlw <- max(logw)
    w <- exp(logw - maxlw)
    # update log likelihood estimate
    ll <- ll + maxlw + log(mean(w))
    normweights <- w / sum(w)
    #
    xhistory[[time+1]] <- xparticles
    whistory[[time+1]] <- normweights
  }
  return(list(ll = ll, xhistory = xhistory, whistory = whistory, ahistory = ahistory))
}