bbricks (version 0.1.1)

HDP2: Create objects of type "HDP2".

Description

Create an object of type "HDP2", which represents the Hierarchical-Dirichlet-Process of two Dirichlet-Process hierarchies: G_m |eta ~ DP(eta,U), m = 1:M G_mj|gamma,G_m ~ DP(gamma,G_m), j = 1:J_m pi_mj|G_mj,alpha ~ DP(alpha,G_mj) z|pi_mj ~ Categorical(pi_mj) k|z,G_mj ~ Categorical(G_mj), if z is a sample from the base measure G_mj u|k,G_m ~ Categorical(G_m), if k is a sample from the base measure G_m theta_u|psi ~ H0(psi) x|theta_u,u ~ F(theta_u) where DP(eta,U) is a Dirichlet Process on positive integers, eta is the "concentration parameter", U is the "base measure" of this Dirichlet process, U is an uniform distribution on all positive integers. DP(gamma,G_m) is a Dirichlet Process on integers with concentration parameter gamma and base measure G_m. DP(alpha,G_mj) is a Dirichlet Process on integers with concentration parameter alpha and base measure G_mj. The choice of F() and H0() can be arbitrary, they are distributions of x and theta_u correspondingly. In the case of HDP2, u, z and k can only be positive integers. This object will be used as a place for recording and accumulating information in the related inference/sampling functions such as posterior(), posteriorDiscard(), MAP() and so on.

Usage

HDP2(
  objCopy = NULL,
  ENV = parent.frame(),
  gamma = list(eta = 1, gamma = 1, alpha = 1, m = 3L, j = c(2, 3, 4), H0aF =
    "GaussianNIW", parH0 = list(m = 0, k = 1, v = 2, S = 1))
)

Arguments

objCopy

an object of type "HDP2". If "objCopy" is not NULL, the function create a new "HDP2" object by copying the content from objCopy, otherwise this new object will be created by using "ENV" and "gamma". Default NULL.

ENV

environment, specify where the object will be created.

gamma

list, a named list of parameters, gamma=list(eta,gamma,alpha,m,j,H0aF,parH0), where gamma$eta is a numeric value specifying the concentration parameter of DP(eta,U), gamma$gamma is a numeric value specifying the concentration parameter of DP(gamma,G_m), gamma$alpha is a numeric value specifying the concentration parameter of DP(alpha,G_mj), gamma$m is the number of groups M, gamma$j is the number of subgroups in each group, must satisfy length(gamma$j)=gamma$m. gamma$H0aF is the name of the "BasicBayesian" object which specifies the structure of H0() and F(). gamma$parH0 is the parameters passed to the selected H0aF. For example, if gamma$H0aF="GaussianNIW", then gamma$parH0 should be a named list of NIW parameters: gamma$parH0=list(m,k,v,S), where gamma$parH0$m is a numeric "location" parameter; gamma$parH0$S is a symmetric positive definite matrix representing the "scale" parameters; gamma$parH0$k and gamma$parH0$v are numeric values.

Value

An object of class "HDP2".

References

Teh, Yee W., et al. "Sharing clusters among related groups: Hierarchical Dirichlet processes." Advances in neural information processing systems. 2005.

See Also

BasicBayesian,GaussianNIW,GaussianNIG,CatDirichlet,CatHDP2,posterior.HDP2,posteriorDiscard.HDP2,marginalLikelihood.HDP2 ...

Examples

# NOT RUN {
## This is an example of Gibbs sampling on a hierarchical mixture model, using HDP2.

## load some hierarchical mixture data, check ?mmhhData for details.
data(mmhhData)
x <- mmhhData$x
ms <- mmhhData$groupLabel
js <- mmhhData$subGroupLabel

## Step1: initialize--------------------------------------------------
maxit <- 50                            #iterative for maxit times
z <- rep(1L,nrow(x))
k <- rep(1L,nrow(x))
u <- rep(1L,nrow(x))
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=2,S=diag(2)*0.001)))
ss <- sufficientStatistics(obj$H,x=x,foreach = TRUE) #sufficient statistics
N <- length(ss)
for(i in 1L:N){                         #initialize z k and u
   tmp <- rPosteriorPredictive(obj = obj,n=1,x=x[i,,drop=FALSE],m=ms[i],j=js[i])
   z[i] <- tmp["z"]
   k[i] <- tmp["k"]
   u[i] <- tmp["u"]
   posterior(obj = obj,ss = ss[[i]],ss1 = u[i],ss2 = k[i],ss3 = z[i],m=ms[i],j = js[i])
}

## Step2: main Gibbs loop---------------------------------------------
it <- 0                                 #iteration tracker
pb <- txtProgressBar(min = 0,max = maxit,style = 3)
while(TRUE){
   for(i in 1L:N){
       ##remove the sample from the posterior info
       posteriorDiscard(obj = obj,ss = ss[[i]],ss1=u[i],ss2=k[i],ss3 = z[i],m=ms[i],j=js[i])
       ##resample a new partition
       tmp <- rPosteriorPredictive(obj = obj,n=1L,x=x[i,,drop=FALSE],m=ms[i],j=js[i])
       z[i] <- tmp["z"]
       k[i] <- tmp["k"]
       u[i] <- tmp["u"]
       posterior(obj = obj,ss = ss[[i]], ss1=u[i],ss2 = k[i],ss3 = z[i],m=ms[i],j=js[i])
   }
   plot(x=x[,1],y=x[,2],col=u)
   it <- it+1
   setTxtProgressBar(pb,it)
   if(it>=maxit){cat("\n");break}
}

# }