## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
#  # install development version from GitHub:
#  # install.packages("devtools")
#  devtools::install_github("chenhaotian/Bayesian-Bricks")

## ----eval=FALSE---------------------------------------------------------------
#  ## Bayesian linear regression
#  
#  library(bbricks)
#  
#  ## lrData is a list of two elements. lrData$x is the sample set of the dependent variable; lrData$X is the sample set of the independent variable
#  ## see ?lrData for details
#  data(lrData)
#  X <- lrData$X                           #a matrix of 1 column
#  x <- lrData$x                           #a numeric vector
#  ## task 1. update the prior into posterior using X and x
#  obj <- GaussianNIG(gamma=list(m=0,V=1,a=1,b=0)) #create a GaussianNIG object
#  ss <- sufficientStatistics(obj = obj,X=X,x=x)   #the sufficient statistics of X and x
#  posterior(obj = obj,ss = ss)                    #add the infomation to the posterior
#  ## task 2. get MAP estimate of beta and sigma^2 from the posterior
#  bsMAP <- MAP(obj)                               #get the MAP estimate of beta and sigma^2
#  bsMAP                                           #print the MAP estimate
#  ## plot the MAP estimate of the regression line
#  plot(X,X%*%bsMAP$betaMAP,type = "l")
#  points(X,x,pch=20)
#  ## task 3. calculate marginal likelihood
#  ## generate some new data
#  Xnew <- matrix(runif(3,min=0,max=),ncol=1)
#  xnew <- Xnew*0.2+rnorm(3,sd=10)
#  marginalLikelihood(obj = obj,X=x,x=x,LOG = TRUE)
#  ## task 4. calculate the posterior prediction
#  ## say we want to predict x at the location X=100
#  predictedSamples <- rPosteriorPredictive(obj = obj,X=matrix(101,ncol = 1),n=1000)
#  ## histogram of the prediction
#  hist(predictedSamples)
#  ## the mean and standard devition of the prediction
#  mean(predictedSamples)
#  sd(predictedSamples)
#  

## ----eval=FALSE---------------------------------------------------------------
#  ## Estimate cancer mortality rates using Gibbs sampling
#  
#  library(bbricks)
#  
#  ## see ?cancerData for details
#  data(cancerData)
#  ## Step1: Initialization----------------------------------------------
#  K <- length(cancerData)                          #then number of cities
#  eta <- 1                                         #assume eta is known, eta=1
#  ## initialize alpha, PI, and sufficient statistics
#  a <- rexp(2,rate = eta)                 #initialize alpha
#  PI <- matrix(0,K,2L)                    #initialize pi
#  cityPrior <- CatDirichlet(gamma = list(alpha=a,uniqueLabels=c("death","no death")))
#  citySS <- lapply(cancerData,function(x){sufficientStatistics(obj = cityPrior,x=x)}) #sufficient statistics of each city
#  ## initialize functions used in Metropolis-Hastings, see ?MetropolisHastings for details
#  ## density of the target distribution
#  dp <- function(a){
#      if(any(a<0)) -Inf
#      else sum(dDir(x=PI,alpha = a,LOG = TRUE))+sum(dexp(x=a,rate = eta,log = TRUE))
#  }
#  ## density of the proposal distribution
#  dq <- function(anew,a){1}                #use a independent proposal
#  ## random sample generator of the proposal distribution
#  rq <- function(x){
#      c(runif(1,x[1]-1,x[1]+1),
#        runif(1,x[2]-1,x[2]+1))
#  }
#  ## Step2: main Gibbs sampling loop between alpha and pi --------------
#  maxit <- 1000
#  burnin <- 500                            #number of burn-in samples
#  meanPI <- numeric(K)                     #place-hoder for the sample mean
#  it <- 1
#  while(it<=maxit){
#      ## Step1: sample pi from p(pi|a,x)-------------
#      for(k in 1L:K){
#          posterior(obj = cityPrior,ss=citySS[[k]])
#          PI[k,] <- rDir(n=1,alpha = cityPrior$gamma$alpha)
#          posteriorDiscard(obj = cityPrior,ss=citySS[[k]])
#      }
#      ## calculate the sample mean
#      if(it>burnin) meanPI <- meanPI+PI[,1]/(maxit-burnin)
#      ## Step2: sample a from p(a|pi,g)--------------
#      ## use Metropolis-Hastings
#      a <- MetropolisHastings(nsamples = 1,xini = a,dp=dp,dq=dq,rq=rq)
#      ## increase iteration counter
#      it <- it+1
#  }
#  ## Step3: plot the result---------------------------------------------
#  ## black bars are the sample mean from the hierarchical Bayesian model
#  ## blue bars are the MLE of the mortality rates.
#  plot(1:K,meanPI,type = "h",xlab = "city",ylab = "mortality rate",lwd=3)
#  lines(1:K+0.2,sapply(cancerData,function(l){sum(l=="death")/length(l)}),type = "h",col = "blue",lwd = 3)
#  legend(1, 0.005, legend=c("Sample Mean", "MLE"),col=c("black", "blue"), lty=c(1,1), cex=1,lwd = 3)

## ----eval=FALSE---------------------------------------------------------------
#  ## Get the MAP estimate of pi and theta using EM algorithm.
#  
#  library(bbricks)
#  
#  ## load some mixture of Gaussian samples.
#  ## mmData is a numeric matrix with 2 columns, each row is a sample
#  ## see ?mmData for details
#  data(mmData)
#  K <- 4L                                 #number of clusters(mixtures components)
#  z <- matrix(runif(nrow(mmData)*K),nrow(mmData),K) #the expected cluster label of each observation
#  allK <- 1L:K    #temp variable, all component labels
#  allZ <- rep(allK,each=nrow(mmData))     #temp variable, all possible cluster labels for all observations
#  ## z, pi and alpha are distributed as a Categorical-Dirichlet sturcture:
#  mc <- CatDirichlet(gamma = list(alpha=0.5,uniqueLabels=allK)) # create a CatDirichlet object to track the posterior info, see ?CatDirichlet for details
#  ## each component distribution is a Gaussian-NIW structure:
#  ec <- replicate(K,GaussianNIW(gamma = list(m=c(0,0),k=0.00001,v=2,S=diag(2)))) # create a GaussianNIW object to track the posterior info of each mixture component, see ?GaussianNIW for details
#  mcMAP <- MAP(mc)                        #initialize the MAP estimate of pi
#  ecMAP <- replicate(K,list(muMAP=runif(2),sigmaMAP=diag(2)),simplify = FALSE) #initialize the MAP estimate of theta
#  ## The main EM loop
#  maxit <- 100                            #number of EM loops
#  it <- 1
#  while(it<=maxit){
#      ## E-step---------------------------------------------------------
#      ## calculate the expected cluster labels: p(z|pi,theta)
#      for(k in allK) z[,k] <- dGaussian(x=mmData,mu = ecMAP[[k]]$muMAP,Sigma=ecMAP[[k]]$sigmaMAP)+log(mcMAP[k])
#      z <- exp(z-logsumexp(z))            #use logsumexp() to avoid numerical underflow
#      ## calculate the expected sufficient statistics
#      ssComponents <- lapply(allK,function(k){
#          sufficientStatistics_Weighted(obj = ec[[k]],x=mmData,w=z[,k])
#      })                                  #the expected sufficient statistics of each Gaussian component
#      ssPi <- sufficientStatistics_Weighted(obj = mc,x=allZ,w=as.vector(z)) #the expected sufficient statistics of the cluster label distribution
#      ## M-step---------------------------------------------------------
#      ## use the sufficient statistics to update the prior distributions:
#      for(k in allK) posterior(obj = ec[[k]],ss=ssComponents[[k]]) #update component distributions
#      posterior(obj = mc,ss = ssPi)                                #update cluster label distribution
#      ## calculate the MAP estimates from posterior:
#      mcMAP <- MAP(mc)
#      ecMAP <- lapply(ec,MAP)
#      ## Reset the priors for next EM loop-----------------------------------------
#      ## to prepare for the next EM iteration, discard the sufficient statistics info from the posteriors:
#      for(k in allK) posteriorDiscard(obj = ec[[k]],ss=ssComponents[[k]])
#      posteriorDiscard(obj = mc,ss = ssPi)
#      ## increase the iteration counter
#      it <- it+1
#  }
#  
#  plot(mmData,col=apply(z,1,which.max)) #plot the best estimates
#  mcMAP                                 #the MAP estimate of pi
#  ecMAP                                 #the MAP estimate of theta_z

## ----eval=FALSE---------------------------------------------------------------
#  ## Sample cluster labels z from DP-MM using Gibbs sampling
#  
#  library(bbricks)
#  
#  ## load some mixture of Gaussian samples.
#  ## mmData is a numeric matrix with 2 columns, each row is a sample
#  ## see ?mmData for details
#  data(mmData)
#  maxit <- 100                            #number of total samples
#  burnin <- 50                            #number of burnin samples
#  ## Step1: Initialization -----------------------------------------
#  obj <- DP(gamma = list(alpha=10,H0aF="GaussianNIW",parH0=list(m=c(0,0),k=0.001,v=2,S=diag(2)))) #create a DP object to track all the changes, the DP object in this case is a combination of a CatDP object and a GaussianNIW object
#  z <- matrix(1L,nrow(mmData),maxit-burnin)    #place-holder for the sampled z
#  ss <- sufficientStatistics(obj,x=mmData,foreach = TRUE) #sufficient statistics of each observed sample
#  N <- nrow(mmData)
#  for(i in 1L:N){ # initialize labels before Gibbs sampling
#      z[i,1] <- rPosteriorPredictive(obj = obj,n=1,x=mmData[i,,drop=FALSE])
#      posterior(obj = obj,ss = ss[[i]], z = z[i,1])
#  }
#  ## Step2: Main Gibbs sampling loop--------------------------------
#  it <- 1                                 #iteration tracker
#  pb <- txtProgressBar(min = 0,max = maxit,style = 3)
#  while(it<=maxit){
#      if(it>burnin) colIdx <- it-burnin
#      else colIdx <- 1
#      for(i in 1L:N){
#          ## remove the sample information from the posterior
#          posteriorDiscard(obj = obj,ss = ss[[i]],z=z[i,colIdx])
#          ## get a new sample
#          z[i,colIdx] <- rPosteriorPredictive(obj = obj,n=1,x=mmData[i,,drop=FALSE])
#          ## add the new sample information to the posterior
#          posterior(obj = obj,ss = ss[[i]],z=z[i,colIdx])
#      }
#      if(it>burnin & colIdx<ncol(z)) z[,colIdx+1] <- z[,colIdx] #copy result of previous iteration
#      it <- it+1
#      setTxtProgressBar(pb,it)
#      if(it>=maxit){cat("\n");break}
#      plot(x=mmData[,1],y=mmData[,2],col=z[,colIdx]) #to see how the labels change in each iteration
#  }
#  ## Step3: Estimate group labels of each observation---------------
#  ## pick the most frequent z as the best estimate
#  zBest <- apply(z,1,function(l){
#      tmp <- table(l)
#      names(tmp)[which.max(tmp)]
#  })
#  plot(x=mmData[,1],y=mmData[,2],col=zBest)

## ----eval=FALSE---------------------------------------------------------------
#  ## 1. add the information of the 4 observed samples to the DP object
#  ssObserved <- sufficientStatistics(obj=obj,x=mmData[c(50,100,150,200),,drop=FALSE],foreach = TRUE)
#  for(i in 1L:4L) posterior(obj = obj,ss = ssObserved[[i]], z = i) # the choice of cluster label 'z' for the 4 observed samples are arbitrary, as long as they are different from each other. In this example I simply use z=1L:4L.
#  ## 2. remove the 4 samples from the upcoming Gibbs sampling procedure
#  mmData <- mmData[-c(50,100,150,200),]

## ----eval=FALSE---------------------------------------------------------------
#  ## Sample cluster labels k from HDP-MM using Gibbs sampling
#  
#  library(bbricks)
#  
#  ## load some mixture of Gaussian samples.
#  ## mmhData is a list of two elements. mmhData$x is a matrix of Gaussian observations, each row is an observation; mmhData$groupLabel is the group label of each observation.
#  ## see ?mmhData for details
#  data(mmhData)
#  x <- mmhData$x
#  js <- mmhData$groupLabel
#  ## Step1: Initialization------------------------------------------
#  maxit <- 50                             #iterative for maxit times
#  burnin <- 30                            #number of burn in samples
#  ## create a HDP object to track all the changes, the HDP object in this case is a combination of a CatHDP object and a GaussianNIW object:
#  obj <- HDP(gamma = list(gamma=1,j=max(js),alpha=1,
#                          H0aF="GaussianNIW",
#                          parH0=list(m=c(0,0),k=0.001,v=2,S=diag(2)*0.01)))
#  ss <- sufficientStatistics(obj$H,x=x,foreach = TRUE) #sufficient statistics
#  set.seed(1)
#  z <- rep(1L,nrow(x))
#  k <- matrix(1L,nrow(x),maxit-burnin)    #place-holder for the sampled k
#  N <- length(ss)
#  for(i in 1L:N){# initialize k and z
#      tmp <- rPosteriorPredictive(obj = obj,n=1,x=x[i,,drop=FALSE],j=js[i])
#      z[i] <- tmp["z"]
#      k[i,1] <- tmp["k"]
#      posterior.HDP(obj = obj,ss = ss[[i]],ss1 = k[i],ss2 = z[i],j = js[i])
#  }
#  ## Step2: main Gibbs loop---------------------------------------------
#  it <- 1                                 #iteration tracker
#  pb <- txtProgressBar(min = 0,max = maxit,style = 3)
#  while(it<=maxit){
#      if(it>burnin) colIdx <- it-burnin
#      else colIdx <- 1
#      for(i in 1L:N){
#          ## remove the sample from the posterior info
#          posteriorDiscard(obj = obj,ss = ss[[i]],ss1=k[i,colIdx],ss2=z[i],j=js[i])
#          ## resample a new partition
#          tmp <- rPosteriorPredictive(obj = obj,n=1,x=x[i,,drop=FALSE],j=js[i])
#          z[i] <- tmp["z"]
#          k[i,colIdx] <- tmp["k"]
#          ## add the information of the new sample
#          posterior(obj = obj,ss = ss[[i]], ss1=k[i,colIdx],ss2 = z[i],j=js[i])
#      }
#      if(it>burnin & colIdx<ncol(k)) k[,colIdx+1] <- k[,colIdx] #copy result of previous iteration
#      it <- it+1
#      plot(x=x[,1],y=x[,2],col=k[,colIdx])         #to visualize the group label dynamics
#      setTxtProgressBar(pb,it)
#  }
#  ## Step3: Estimate group labels of each observation---------------
#  ## pick the most frequent k as the best estimate
#  kBest <- apply(k,1,function(l){
#      tmp <- table(l)
#      names(tmp)[which.max(tmp)]
#  })
#  plot(x=x[,1],y=x[,2],col=kBest)

## ----eval=FALSE---------------------------------------------------------------
#  ## HDP-LDA on the farm-ads corpus
#  
#  ## load a subset of farm ads data from https://archive.ics.uci.edu/ml/datasets/Farm+Ads
#  ## see ?farmadsData for details
#  data(farmadsData)
#  word <- farmadsData$word
#  document <- farmadsData$document
#  ## Step1: Initialization------------------------------------------
#  set.seed(1)
#  maxit <- 30                            #iterative for maxit times
#  z <- rep(1L,length(word))
#  k <- rep(1L,length(word))
#  ## initialize
#  uniqueWords <- unique(word)
#  obj <- HDP(gamma = list(gamma=1,j=max(document),alpha=1,H0aF="CatDirichlet",parH0=list(alpha=rep(0.5,length(uniqueWords)),uniqueLabels=uniqueWords))) #create a HDP object to track all the changes, the HDP object in this case is a combination of a CatHDP object and a CatDrirchlet object
#  N <- length(word)
#  ## initialize k and z
#  for(i in 1L:N){
#      tmp <- rPosteriorPredictive(obj = obj,n=1,x=word[i],j=document[i])
#      z[i] <- tmp["z"]
#      k[i] <- tmp["k"]
#      posterior(obj = obj,ss = word[i], ss2 = z[i],j=document[i],ss1=k[i])
#  }
#  ## Step2: main Gibbs loop---------------------------------------------
#  it <- 1                                 #iteration tracker
#  pb <- txtProgressBar(min = 0,max = maxit,style = 3)
#  while(it<=maxit){
#      for(i in 1L:N){
#          posteriorDiscard.HDP(obj = obj,ss = word[i],ss1=k[i],ss2=z[i],j=document[i]) #remove the sample information from the posterior
#          tmp <- rPosteriorPredictive(obj = obj,n=1,x=word[i],j=document[i])   #get a new sample
#          z[i] <- tmp["z"]
#          k[i] <- tmp["k"]
#          posterior(obj = obj,ss = word[i],ss1=k[i], ss2 = z[i],j=document[i]) #add the new sample information to the posterior
#      }
#      it <- it+1
#      setTxtProgressBar(pb,it)
#  }
#  ## Step3: plot the result --------------------------------------------
#  ## see which topics are most frequently appeared:
#  order(sapply(obj$X,function(l){sum(l$gamma$alpha)}),decreasing = TRUE)
#  ## seems topic 2 and 1 appear the most, let's plot them:
#  ## install.packages("wordcloud") # for word-cloud
#  ## install.packages("RColorBrewer") # color palettes
#  ## print topic 1
#  wordcloud:: wordcloud(words = obj$X[[1]]$gamma$uniqueLabels,
#                        freq = obj$X[[1]]$gamma$alpha,
#                        min.freq = 1,
#                        max.words=100,
#                        random.order=FALSE, rot.per=0.35,
#                        colors=RColorBrewer::brewer.pal(5, "Set1"))
#  ## print topic 2
#  wordcloud:: wordcloud(words = obj$X[[2]]$gamma$uniqueLabels,
#                        freq = obj$X[[2]]$gamma$alpha,
#                        min.freq = 1,
#                        max.words=100,
#                        random.order=FALSE, rot.per=0.35,
#                        colors=RColorBrewer::brewer.pal(5, "Set1"))

## ----eval=FALSE---------------------------------------------------------------
#  ## Sample cluster labels u from HDP2-MM using Gibbs sampling
#  
#  library(bbricks)
#  
#  ## load some mixture of Gaussian samples.
#  ## mmhData is a list of 3 elements. mmhhData$x is a matrix of Gaussian observations, each row is an observation; mmhhData$groupLabel is the group label of each observation. mmhhData$subGroupLabel is the subgroup label of each observation.
#  ## there are 2 groups, group1 has 10 subgroups, group has 20 subgroups.
#  ## see ?mmhhData for details
#  data(mmhhData)
#  x <- mmhhData$x
#  groupLabel <- mmhhData$groupLabel
#  subGroupLabel <- mmhhData$subGroupLabel
#  ## Step1: Initialization------------------------------------------
#  maxit <- 50                            #iterative for maxit times
#  burnin <- 20                           #number of burn in samples
#  ## create a HDP2 object to track all the changes, the HDP2 object in this case is a combination of a CatHDP2 object and a GaussianNIW object
#  obj <- HDP2(gamma = list(eta=1,gamma=1,alpha=1,m=2L,j=c(10L,20L),H0aF="GaussianNIW",parH0=list(m=c(0,0),k=0.001,v=1.1,S=diag(2)*0.001)))
#  ss <- sufficientStatistics(obj$H,x=x,foreach = TRUE) #sufficient statistics
#  z <- rep(1L,nrow(x))
#  k <- rep(1L,nrow(x))
#  u <- matrix(1L,nrow(x),maxit-burnin)
#  N <- length(ss)
#  ## initialization k, z and u
#  for(i in 1L:N){
#      tmp <- rPosteriorPredictive(obj = obj,n=1,x=x[i,,drop=FALSE],m=groupLabel[i],j=subGroupLabel[i])
#      z[i] <- tmp["z"]
#      k[i] <- tmp["k"]
#      u[i,1] <- tmp["u"]
#      posterior.HDP2(obj = obj,ss = ss[[i]],ss1 = u[i,1],ss2 = k[i],ss3 = z[i],m=groupLabel[i],j = subGroupLabel[i])
#  }
#  ## Step2: main Gibbs loop---------------------------------------------
#  it <- 1                                 #iteration counter
#  pb <- txtProgressBar(min = 0,max = maxit,style = 3)
#  while(it<=maxit){
#      if(it>burnin) colIdx <- it-burnin
#      else colIdx <- 1
#      for(i in 1L:N){
#          ## remove the sample from the posterior info
#          posteriorDiscard(obj = obj,ss = ss[[i]],ss1=u[i,colIdx],ss2=k[i],ss3 = z[i],m=groupLabel[i],j=subGroupLabel[i])
#          ## resample a new partition
#          tmp <- rPosteriorPredictive(obj = obj,n=1L,x=x[i,,drop=FALSE],m=groupLabel[i],j=subGroupLabel[i])
#          z[i] <- tmp["z"]
#          k[i] <- tmp["k"]
#          u[i,colIdx] <- tmp["u"]
#          ## add the information of the new sample
#          posterior(obj = obj,ss = ss[[i]], ss1=u[i,colIdx],ss2 = k[i],ss3 = z[i],m=groupLabel[i],j=subGroupLabel[i])
#      }
#      if(it>burnin & colIdx<ncol(u)) u[,colIdx+1] <- u[,colIdx] #copy result of previous iteration
#      it <- it+1
#      plot(x=x[,1],y=x[,2],col=u[,colIdx])
#      setTxtProgressBar(pb,it)
#  }
#  ## Step3: Estimate group labels of each observation---------------
#  ## pick the most frequent u as the best estimate
#  uBest <- apply(u,1,function(l){
#      tmp <- table(l)
#      names(tmp)[which.max(tmp)]
#  })
#  plot(x=x[,1],y=x[,2],col=uBest)

