Description Usage Arguments Value References See Also Examples
View source: R/Dirichlet_Process.r
Create an object of type "HDP2", which represents the Hierarchical-Dirichlet-Process with two Dirichlet-Process hierarchies:
G |eta \sim DP(eta,U)
G_m|gamma,G \sim DP(gamma,G), m = 1:M
pi_{mj}|G_m,alpha \sim DP(alpha,G_m), j = 1:J_m
z|pi_{mj} \sim Categorical(pi_{mj})
k|z,G_m \sim Categorical(G_m),\textrm{ if z is a sample from the base measure }G_m
u|k,G \sim Categorical(G),\textrm{ if k is a sample from the base measure G}
theta_u|psi \sim H0(psi)
x|theta_u,u \sim 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) is a Dirichlet Process on integers with concentration parameter gamma and base measure G. DP(alpha,G_m) is a Dirichlet Process on integers with concentration parameter alpha and base measure G_m. The choice of F() and H0() can be described by an arbitrary "BasicBayesian" object such as "GaussianGaussian","GaussianInvWishart","GaussianNIW", "GaussianNIG", "CatDirichlet", and "CatDP". See ?BasicBayesian
for definition of "BasicBayesian" objects, and see for example ?GaussianGaussian
for specific "BasicBayesian" instances. As a summary, An "HDP2" object is simply a combination of a "CatHDP2" object (see ?CatHDP2
) and an object of any "BasicBayesian" type.
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.
1 2 3 4 5 6 |
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), gamma$alpha is a numeric value specifying the concentration parameter of DP(alpha,G_m), 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. |
An object of class "HDP2".
Teh, Yee W., et al. "Sharing clusters among related groups: Hierarchical Dirichlet processes." Advances in neural information processing systems. 2005.
BasicBayesian
,GaussianNIW
,GaussianNIG
,CatDirichlet
,CatHDP2
,posterior.HDP2
,posteriorDiscard.HDP2
,marginalLikelihood.HDP2
...
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 | ## 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}
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.