R/DPcluster.R

In morner75/BNP:

#' Conditional distribution of cluster membership indicators
#'
#' @param s an integer vector / a sequence of cluster memebership indicators
#' @param M a numeric / precision parameter
#' @return a list of one numeric vector and one scalar / conditional dist. and probability of p(s_i | s_{-1})
# debugging:
# S <-SS[1:4] #c(5,3,3,2,3,4,5,4,5,6,4,4,3,3,7) ; M=1
cdclust <- function(S, M=1){
  n <- length(S)
  S <- as.factor(S)
  clust.val <- unique(S[-n])
  dist <- c(vapply(clust.val, function(i) sum(S[-n]==i)/(length(S)+M-1),FUN.VALUE = numeric(1)), M/(M+length(S)-1))
  names(dist) <- c(format(clust.val,digits = 3),"new")
  if(S[n]%in%clust.val){
    prob <- dist[c(clust.val==S[n], FALSE)]
  }else{
    prob <- dist[length(dist)]
  }
  list(dist=dist,prob=prob)
}


#' Marginal distribution of cluster membership indicators
#'
#' @param S an integer vector / a sequence of cluster memebership indicators
#' @param M a numeric / precision parameter
#' @return alist of one numeric vector and one scalar / conditional dist. and probability of p(s_i | s_{-1})
# debugging:
# S <-c(5,3,3,2,3,4,5,4,5,6,4,4,3,3,7) ; M=1 ; log = TRUE
dclust <- function(S, M=1, log = FALSE){
  prob <- prod(vapply(seq_along(S),FUN=function(i) cdclust(S[1:i],M=M)$prob, FUN.VALUE=numeric(1)))
  if(log==TRUE) prob <- log(prob)
  prob
}