R/approximatesmoothC++.R

In FVDDPpkg: Implement Fleming-Viot-Dependent Dirichlet Processes

#' Approximate the smoothing distribution of a Fleming-Viot latent signal
#'
#' @param fvddp.past An instance of class `fvddp`, progressively updated ad propagated
#' with data referring to past times via [FVDDPpkg::update()] and [FVDDPpkg::propagate()]
#' (or its approximate version, [FVDDPpkg::approx.propagate()]).
#' @param fvddp.future Same as `fvddp.past`, but in this case the propagation has
#' been performed with time data from times later than the one to be estimated. Its
#' hyperparameters must be equals to the ones of `fvddp.past`.
#' @param t.past The time between the last collection of data (in the past) and the
#' time at which the smoothing is performed.
#' @param t.future Same as `t.past`, but in this case it is referred to the future.
#' `t.future` is positive too.
#' @param y.new The data collected at the time the smoothing is performed.
#' @param N the amount of samples to be drawn in order to perform the approximation.
#'
#' @return An object of class `fvddp`, with the same hyperparmeters as `fvddp.past`
#' and `fvddp.future`. Since this function is a Monte-Carlo based
#' approximation of [FVDDPpkg::smooth()], the outputs are similar.
#' @export
#'
#' @examples
#' FVDDP = initialize(3, function(x) rbinom(x, 10, .2),
#'                    function(x) dbinom(x, 10, .2), TRUE)
#' FVDDP.PAST = update(FVDDP, c(2,3))
#' FVDDP.FUTURE = update(FVDDP, c(4))
#' FVDDP.FUTURE = propagate(FVDDP.FUTURE, 0.5)
#' FVDDP.FUTURE = update(FVDDP.FUTURE, c(1))
#' approx.smooth(fvddp.past = FVDDP.PAST, fvddp.future = FVDDP.FUTURE,
#'               t.past = 0.4, t.future = 0.3, y.new = c(1,3), N = 20000)
#'
#' @seealso [smooth()] for a (slower) exact computation
approx.smooth = function(fvddp.past, fvddp.future, t.past, t.future, y.new, N){

  #check the class of the arguments
  if (!inherits(fvddp.past, "fvddp")) stop(deparse(substitute(fvddp.past)),
                                         ' not in "fvddp" class')
  #check the class of the fvddp
  if (!inherits(fvddp.future, "fvddp")) stop(deparse(substitute(fvddp.future)),
                                           ' not in "fvddp" class')

  #we need the processes (from both past and future) to have the same theta and P0
  if ((fvddp.past$theta != fvddp.future$theta) |
      (fvddp.past$is.atomic != fvddp.past$is.atomic)) stop('Unmatched parameters')

  #initialize a brand new fvddp
  else {
    fvddp = initialize(theta = fvddp.past$theta, sampling.f = fvddp.past$P0.sample,
                       density.f = fvddp.past$P0.density, atomic = fvddp.past$is.atomic)
  }

  #rename the matrices of multiplicities
  M.past = fvddp.past$M
  M.future = fvddp.future$M

  #get the final value of y.star
  y.star = sort(unique(c(fvddp.past$y.star, fvddp.future$y.star, unique(y.new))))

  #check if some values need to be added on y.star on the past or on the future
  new.y.star.past = setdiff(y.star, fvddp.past$y.star)
  new.y.star.future = setdiff(y.star, fvddp.future$y.star)

  #if it is necessary, add some columns of zeros on M.past
  if (length(new.y.star.past) != 0) {
    M.past = cbind(M.past, matrix(0, nrow=nrow(M.past), ncol=length(new.y.star.past),
                                  dimnames = list(NULL,new.y.star.past)))
    M.past = M.past[,order(as.numeric(colnames(M.past)))]

    #modify the type, if it has been coverted to vector
    if (is.null(dim(M.past))) M.past = matrix(M.past, ncol=length(y.star))
  }

  #if it is necessary, add some columns of zeros on M.past
  if (length(new.y.star.future) !=  0){
    M.future = cbind(M.future, matrix(0, nrow=nrow(M.future), ncol=length(new.y.star.future),
                                      dimnames = list(NULL,new.y.star.future)))
    M.future = M.future[,order(as.numeric(colnames(M.future)))]

    #modify the type, if it has been coverted to vector
    if (is.null(dim(M.future))) M.future = matrix(M.future, ncol=length(y.star))
  }

  #write n as a vector of multiplicities
  n  = table(factor(y.new, levels=y.star))

  #save the vectors of weights
  w.past = fvddp.past$w
  w.future = fvddp.future$w

  #keep these useful informations
  n.past.max = max(rowSums(M.past))
  n.future.max = max(rowSums(M.future))
  tot = max(n.past.max, n.future.max)
  theta.P0.star = fvddp$theta*fvddp$P0.density(y.star)


  #these objects are necessary to keep in memory some terms
  lambda = (0:tot)*(fvddp$theta + (0:tot) -1)/2
  nonatomic.w.mat = matrix(nrow = n.past.max + 1, ncol = n.future.max +1)


  lst = montecarlo_sample_smooth_cpp(M.past, M.future, n, t.past, t.future, w.past,
                                     w.future, N, fvddp$is.atomic, lambda, fvddp$theta,
                                     theta.P0.star, nonatomic.w.mat)

  fvddp$y.star = y.star
  fvddp$M = lst$M
  colnames(fvddp$M) = fvddp$y.star
  fvddp$w = lst$w/sum(lst$w)

  #assign the fvddp to the appropriate class
  class(fvddp) = 'fvddp'

  #return the fvddp list
  return(fvddp)
}