R/posteriorsampling.R

In FVDDPpkg: Implement Fleming-Viot-Dependent Dirichlet Processes

#' Draw values from the posterior distribution FVDDP
#'
#' @param fvddp The instance of class `fvddp` the values are drawn from.
#' @param N The amount of values to draw.
#'
#' @return A vector of length `N` of values drawn either from the centering of the
#' FVDDP (the input) or from the empirical probability measure generated by past
#' observations. The difference between this function and [FVDDPpkg::predictive.struct()]
#' is that in this case the process is not update with respect to any drawn value.
#' @export
#'
#' @examples
#' #create a dummy process and sample some values from it
#' FVDDP = initialize(7, function(x) rbeta(x, 3,3),
#'                    function(x) dgamma(x, 3,3), FALSE)
#' FVDDP = update(FVDDP, rep(0:1, 2))
#' posterior.sample(fvddp = FVDDP, N = 100)
posterior.sample = function(fvddp, N) {

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

  #initialize the empty vector of new values
  y = rep(NA, N)

  #rename these object for practical purpose
  M = fvddp$M
  w = fvddp$w

  #compute the probabilities q for each y_j
  q = t(M * matrix(rep(1/(rowSums(M)+fvddp$theta), ncol(M)),ncol=ncol(M))) %*% fvddp$w

  #q[K+1] whas what we called q_0
  K = length(q)
  q[K+1] = 1 - sum(q)

  #get the indexes, and find those diffent from q_0
  j = sample.int(K+1, size=N, replace=T, prob=q)
  idx = j != (K+1)

  #subsitute the values of y.star in y
  y[idx] = fvddp$y.star[j[idx]]

  #draw the new values from P0
  y[!idx] = fvddp$P0.sample(sum(!idx))

  #return the value
  return(y)
}