rDPGibbs               package:bayesm               R Documentation

_D_e_n_s_i_t_y _E_s_t_i_m_a_t_i_o_n _w_i_t_h _D_i_r_i_c_h_l_e_t _P_r_o_c_e_s_s _P_r_i_o_r _a_n_d _N_o_r_m_a_l _B_a_s_e

_D_e_s_c_r_i_p_t_i_o_n:

     'rDPGibbs' implements a Gibbs Sampler to draw from the posterior
     for a normal mixture problem with a Dirichlet Process prior.  A
     natural conjugate base prior is used along with priors on the
     hyper  parameters of this distribution. One interpretation of this
     model is as a normal mixture with a random number of components
     that can grow with the sample size.

_U_s_a_g_e:

     rDPGibbs(Prior, Data, Mcmc)

_A_r_g_u_m_e_n_t_s:

   Prior: list(Prioralpha,lambda_hyper) 

    Data: list(y) 

    Mcmc: list(R,keep,maxuniq,SCALE,gridsize) 

_D_e_t_a_i_l_s:

     Model: 
      y_i ~ N(mu_i,Sigma_i). 

     Priors:
      theta_i=(mu_i,Sigma_i) ~ DP(G_0(lambda),alpha)
      G_0(lambda):
      mu_i | Sigma_i ~ N(0,Sigma_i (x) a^{-1})
      Sigma_i ~ IW(nu,nu*v*I)

     lambda(a,nu,v):
      a ~ uniform on grid[alim[1],alimb[2]]
      nu ~ uniform on grid[dim(data)-1 + exp(nulim[1]),dim(data)-1
     +exp(nulim[2])]
      v ~ uniform on grid[vlim[1],vlim[2]]

     alpha ~ (1-(alpha-alphamin)/(alphamax-alphamin))^power 
      alpha= alphamin then expected number of components = Istarmin 
      alpha= alphamax then expected number of components = Istarmax 

     list arguments

     Data:

   '_y' N x k matrix of observations on k dimensional data

     Prioralpha:

   '_I_s_t_a_r_m_i_n' expected number of components at lower bound of support
        of alpha

   '_I_s_t_a_r_m_a_x' expected number of components at upper bound of support
        of alpha

   '_p_o_w_e_r' power parameter for alpha prior

     lambda_hyper:

   '_a_l_i_m' defines support of a distribution,def:c(.01,2) 

   '_n_u_l_i_m' defines support of nu distribution, def:c(.01,3) 

   '_v_l_i_m' defines support of v distribution, def:c(.1,4) 

     Mcmc:

   '_R' number of mcmc draws

   '_k_e_e_p' thinning parm, keep every keepth draw

   '_m_a_x_u_n_i_q' storage constraint on the number of unique components

   '_S_C_A_L_E' should data be scaled by mean,std deviation before posterior
        draws, def: TRUE

   '_g_r_i_d_s_i_z_e' number of discrete points for hyperparameter priors,def:
        20

     output:
      the basic output are draws from the predictive distribution of
     the data in the object, 'nmix'.  The average of these draws is the
     Bayesian analogue of a density estimate.

     nmix:

   '_p_r_o_b_d_r_a_w' R/keep x 1 matrix of 1s

   '_z_d_r_a_w' R/keep x N matrix of draws of indicators of which component
        each obs is assigned to

   '_c_o_m_p_d_r_a_w' R/keep list of draws of normals

     Output of the components is in the form of a list of lists. 
      compdraw[[i]] is ith draw - list of lists. 
      compdraw[[i]][[1]] is list of parms for a draw from predictive. 
      compdraw[[i]][1]][[1]] is the mean vector.
     compdraw[[i]][[1]][[2]] is the inverse of Cholesky root. Sigma =
     t(R)%*%R, R^{-1} = compdraw[[i]][[1]][[2]].

_V_a_l_u_e:

    nmix: a list containing: probdraw,zdraw,compdraw

alphadraw: vector of draws of DP process tightness parameter

  nudraw: vector of draws of base prior hyperparameter

   adraw: vector of draws of base prior hyperparameter

   vdraw: vector of draws of base prior hyperparameter

_N_o_t_e:

     we parameterize the prior on Sigma_i such that mode(Sigma)=
     nu/(nu+2) vI. The support of nu enforces valid IW density;
     nulim[1] > 0

     We use the structure for 'nmix' that is compatible with the
     'bayesm' routines for finite mixtures of normals. This allows us
     to use the same summary and plotting methods.  

     The default choices of alim,nulim, and vlim determine the location
     and approximate size of candidate "atoms" or possible normal
     components. The defaults are sensible given that we scale the
     data.  Without  scaling, you want to insure that alim is set for a
     wide enough range of values (remember a is a precision parameter)
     and the v is big enough to propose Sigma matrices wide enough to
     cover the data range.  

     A careful analyst should look at the posterior distribution of a,
     nu, v to make sure that the support is set correctly in alim,
     nulim, vlim.  In other words, if we see the posterior bunched up
     at one end of these support ranges, we should widen the range and
     rerun.  

     If you want to force the procedure to use many small atoms, then
     set nulim to consider only large values and  set vlim to consider
     only small scaling constants.  Set Istarmax to a large number. 
     This will create a very "lumpy" density estimate somewhat like the
     classical Kernel density estimates. Of course, this is not advised
      if you have a prior belief that densities are relatively smooth.

_A_u_t_h_o_r(_s):

     Peter Rossi, Graduate School of Business, University of Chicago,
     Peter.Rossi@ChicagoGsb.edu.

_S_e_e _A_l_s_o:

     'rnmixGibbs','rmixture', 'rmixGibbs' , 'eMixMargDen', 'momMix',
     'mixDen', 'mixDenBi'

_E_x_a_m_p_l_e_s:

     if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

     ## simulate univariate data from Chi-Sq

     set.seed(66)
     N=200
     chisqdf=8; y1=as.matrix(rchisq(N,df=chisqdf))

     ## set arguments for rDPGibbs

     Data1=list(y=y1)
     Prioralpha=list(Istarmin=1,Istarmax=10,power=.8)
     Prior1=list(Prioralpha=Prioralpha)

     Mcmc=list(R=R,keep=1,maxuniq=200)

     out1=rDPGibbs(Prior=Prior1,Data=Data1,Mcmc)

     if(0){
     ## plotting examples
     rgi=c(0,20); grid=matrix(seq(from=rgi[1],to=rgi[2],length.out=50),ncol=1)
     deltax=(rgi[2]-rgi[1])/nrow(grid)
     plot(out1$nmix,Grid=grid,Data=y1)
     ## plot true density with historgram
     plot(range(grid[,1]),1.5*range(dchisq(grid[,1],df=chisqdf)),type="n",xlab=paste("Chisq ; ",N," obs",sep=""), ylab="")
     hist(y1,xlim=rgi,freq=FALSE,col="yellow",breaks=20,add=TRUE)
     lines(grid[,1],dchisq(grid[,1],df=chisqdf)/(sum(dchisq(grid[,1],df=chisqdf))*deltax),col="blue",lwd=2)
     }

     ## simulate bivariate data from the  "Banana" distribution (Meng and Barnard) 
     banana=function(A,B,C1,C2,N,keep=10,init=10)
     { R=init*keep+N*keep
        x1=x2=0
        bimat=matrix(double(2*N),ncol=2)
       for (r in 1:R)
       { x1=rnorm(1,mean=(B*x2+C1)/(A*(x2^2)+1),sd=sqrt(1/(A*(x2^2)+1)))
       x2=rnorm(1,mean=(B*x2+C2)/(A*(x1^2)+1),sd=sqrt(1/(A*(x1^2)+1)))
       if (r>init*keep && r%%keep==0) {mkeep=r/keep; bimat[mkeep-init,]=c(x1,x2)} }
     return(bimat)
     }

     set.seed(66)
     nvar2=2
     A=0.5; B=0; C1=C2=3
     y2=banana(A=A,B=B,C1=C1,C2=C2,1000)

     Data2=list(y=y2)
     Prioralpha=list(Istarmin=1,Istarmax=10,power=.8)
     Prior2=list(Prioralpha=Prioralpha)
     Mcmc=list(R=R,keep=1,maxuniq=200)

     out2=rDPGibbs(Prior=Prior2,Data=Data2,Mcmc)

     if(0){
     ## plotting examples

     rx1=range(y2[,1]); rx2=range(y2[,2])
     x1=seq(from=rx1[1],to=rx1[2],length.out=50)
     x2=seq(from=rx2[1],to=rx2[2],length.out=50)
     grid=cbind(x1,x2)

     plot(out2$nmix,Grid=grid,Data=y2)

     ## plot true bivariate density
     tden=matrix(double(50*50),ncol=50)
     for (i in 1:50){ for (j in 1:50) 
           {tden[i,j]=exp(-0.5*(A*(x1[i]^2)*(x2[j]^2)+(x1[i]^2)+(x2[j]^2)-2*B*x1[i]*x2[j]-2*C1*x1[i]-2*C2*x2[j]))}
     }
     tden=tden/sum(tden)
     image(x1,x2,tden,col=terrain.colors(100),xlab="",ylab="")
     contour(x1,x2,tden,add=TRUE,drawlabels=FALSE)
     title("True Density")
     }

