R/postalphak.R

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

#' Full conditional for alphak
#'
#' @param clusteri A scalar with the index of a cluster
#' @param Y A matrix M x N with the data sequences
#' @param M A scalar representing the number of points available for each data sequence
#' @param K A vector containing the number of change points for each cluster (or its initial values)
#' @param Tl A list containing a vector for each cluster determining the change-point positions in each cluster
#'    (or its initial values)
#' @param cluster A vector containing the cluster assignments for the data sequences (or its initial values)
#' @param sigma2 A vector with the variances of the data sequences (or its initial values)
#'
#' @returns A numerical vector of size `K` + 1 with sampled values from the full conditional of alphak for a given cluster `clusteri`
#'
#' @note
#' This function is called within the Gibbs sampler, but it can be called separately as well.
#'
#' @seealso [gibbs_alg()]
#'
#' @examples
#' data(data)
#' postalphak(M = 50, Y = data, sigma2 = 0.05, K = c(0, 0), Tl = c(50, 50),
#'  cluster = c(1,1,2,1,2), clusteri = 1)
#'
#' @export
#'
postalphak <- function(M, Y, sigma2, K, Tl, cluster, clusteri){
  # filtering the cells in cluster i
  cellsn <- which(cluster == clusteri)

  # selecting the Y values for the cells in cluster i
  Yn <- as.matrix(Y[,cellsn])

  # selecting sigma2 for cells in cluster i
  sigma2n <- sigma2[cellsn]

  # number of change points for cluster i
  k <- K[clusteri]

  # Change point positions for cluster i
  Tln <- Tl[[clusteri]]

  Tlnn <- c()
  if(k == 0){
    Tlnn <- M
    Xn <- as.matrix(rep(1, Tlnn))
  } else{
    Xn <- apply(diag(k+1),2,function(x){rep(x,diff(Tln))})
  }

  Sigman <- solve((sum(sigma2n^(-1)))*t(Xn)%*%Xn)

  if(k == 0){
    mun <- Sigman*(sum(sapply(1:length(cellsn), function(x){(sigma2n[x]^(-1))*t(Xn)%*%Yn[,x]})))
    alphak <- stats::rnorm(1,mun, sqrt(Sigman))
  } else{
    mun <- Sigman%*%as.matrix(rowSums(sapply(1:length(cellsn), function(x){(sigma2n[x]^(-1))*t(Xn)%*%Yn[,x]})))
    alphak <- sapply(1:(k+1), function(x){stats::rnorm(1, mun[x,1], sqrt(diag(Sigman)[x]))})
  }
  return(alphak)
}