R/postK_mk.R

In BayesCPclust: A Bayesian Approach for Clustering Constant-Wise Change-Point Data

#' Marginal probability of K per bin
#'
#' @param k A scalar for the number of changes points in a cluster
#' @param mk A matrix with all possible values to distribute between change points
#' @param m0 A scalar for the number of positions available to define change-points positions
#' @param Yn A vector or matrix with data sequences for a cluster
#' @param sigma2n A vector with the variance of the data sequences in a cluster
#' @param cellsn A vector with the indices of the data sequences in a cluster
#' @param Cr A scalar with the number of data sequences in a cluster
#' @param M A scalar representing the number of points available for each data sequence
#' @param w A scalar representing the minimum number of points in each interval between two change points
#'
#' @returns `postK_mk` returns a numerical value representing the non-normalized probability for a given bin, given k, and a given cluster
#'
#' @note
#' This function is called within [postK()]. It should not be called alone.
#'
#' @seealso [postK()], [gibbs_alg()]
#' @examples
#' data(data)
#' M <- 50; k <- 0; w <- 10;
#' m0 <- M - 1 -(k+1)*w
#' for(k in 0:2){
#' mk <- RcppAlgos::permuteGeneral(0:m0, k + 1,
#' constraintFun = "sum",
#' comparisonFun = "==", limitConstraints = m0,
#' repetition = TRUE)}
#' out <- postK_mk(k = 0, m0 = m0, w = 10, M = 50, Yn = data[,c(1,2,4)],
#'  sigma2n = rep(0.05, 3), cellsn = c(1,2,4), mk = mk[1,], Cr = 3)
#' @export
#'
postK_mk <- function(k, m0, w, M, Yn, sigma2n, cellsn, mk, Cr){
  Tl <- c()
  if(k == 0){
    Tl <- M
    Xn <- as.matrix(rep(1,Tl))

    Vn <- t(Xn)%*%Xn

    A <- sum(sapply(1:length(cellsn), function(x){(sigma2n[x]^(-1))*t(Xn)%*%Yn[,x]}))

  }else{
    for(i in 2:(k+1)){
      Tl[1] <- 1
      Tl[i] <- mk[i-1] + Tl[i-1] + w
    }
    Tl <- append(Tl, M + 1)
    Xn <- apply(diag(k + 1), 2, function(x){rep(x,diff(Tl))})


    Vn <- t(Xn)%*%Xn

    A <- as.matrix(rowSums(sapply(1:length(cellsn), function(x){(sigma2n[x]^(-1))*t(Xn)%*%Yn[,x]})))

  }



  H <- -(M/2)*sum(log(sigma2n)) - ((M*Cr - k - 1)/2)*log(2*pi) - ((k+1)/2)*log(sum(sigma2n^(-1))) - (1/2)*log(det(Vn)) - (1/2)*sum(sapply(1:length(cellsn), function(x){(sigma2n[x]^(-1))*t(Yn[,x])%*%Yn[,x]})) + (1/2)*t(A)%*%solve((sum(sigma2n^(-1)))*t(Xn)%*%Xn)%*%A



  return(H)
}