R/DPM.HODC.R

Defines functions DPM.HODC

Documented in DPM.HODC

#'Hierachical ordered density clustering (HODC) Algorithm with input generated by DPdensity
#'@param pvalue a vector of p-values obtained from large scale statistical hypothesis testing
#'@param v number of iterations set for DPM fitting by "DPdensity"
#'@param DPM.mcmc list
#'@param DPM.prior list
#'@details Without the information of networking, we can have an approximation to the marginal density by DPM model fitting on \strong{r}. Suppose the number of finite mixture normals is equal to L_0+L_1, which means the number of classes we have, we apply HODC algorithm in partitioning the $L_0$ and $L_1$ components into two classes.
#'For this function, the input is generated by Mclust
#'@return a list of HODC algorithm returned parameters.
#'\describe{
#'\item{mean}{the means of each of two cluster for every DPM fitting by "DPdensity"
#'\describe{
#'\item{mu0}{the means of the cluster with smaller mean}
#'\item{mu1}{the means of the cluster with larger mean}
#'}
#'}
#'\item{variance}{the variance of each of two cluster for every DPM fitting by "DPdensity"
#'\describe{
#'\item{var0}{the variances of the cluster with smaller mean}
#'\item{var1}{the variances of the cluster with larger mean}
#'}
#'}
#'\item{probability}{the probability of each of two cluster for every DPM fitting by "DPdensity"
#'\describe{
#'\item{pro0}{the probabilities of the cluster with smaller mean}
#'\item{pro1}{the probabilities of the cluster with larger mean}
#'}
#'}
#'\item{classification}{The classification corresponding to each cluster for every DPM fitting by "DPdensity"}
#'}
#'@examples 
#'\dontrun{
#'###random make the density
#'rstat=c(rnorm(50,mean=1),rnorm(50,mean=2),rnorm(100,mean=4),rnorm(100,mean=8))
#'###transformed into pvalue
#'pvalue=pnorm(-rstat)
#'dpdensityHODC=DPM.HODC(v=5,pvalue)
#'}
#'@export
DPM.HODC<-function(v,pvalue,DPM.mcmc=list(nburn=2000,nsave=1,nskip=0,ndisplay=10),DPM.prior=list(a0=2,b0=1,m2=rep(0,1),s2=diag(100000,1),
                                                                                                 psiinv2=solve(diag(0.5,1)),
                                                                                                 nu1=4,nu2=4,tau1=1,tau2=100)){
  results<-list()
  rstat=Transfer(pvalue)
  
  #nburn <-2000
  #nsave <-1
  #nskip <-0
  #ndisplay <-10
  #mcmc <- list(nburn=nburn,nsave=nsave,nskip=nskip,ndisplay=ndisplay)
  mcmc<-DPM.mcmc
  #prior1 <- list(alpha=0.5,m1=0.2,nu1=100,psiinv1=0.1,k0=100)
  #prior2 <- list(a0=2,b0=1,m2=rep(0,1),s2=diag(100000,1),
  #              psiinv2=solve(diag(0.5,1)),
  #             nu1=4,nu2=4,tau1=1,tau2=100)
  
  fit<-DPpackage::DPdensity(rstat,prior=DPM.prior,mcmc=mcmc,status=TRUE)
  
  for(oo in 1:v){
    cat("iter: ", oo,"\n")
    utils::flush.console()
    nburn <-0
    nsave <-1
    nskip <-0
    ndisplay <-10
    mcmc.conti <- list(nburn=nburn,nsave=nsave,nskip=nskip,ndisplay=ndisplay)
    #prior1 <- list(alpha=0.5,m1=0.2,nu1=100,psiinv1=0.1,k0=100)
    #prior2 <- list(a0=2,b0=1,m2=rep(0,1),s2=diag(100000,1),
    #              psiinv2=solve(diag(0.5,1)),
    #             nu1=4,nu2=4,tau1=1,tau2=100)
    fit<-DPpackage::DPdensity(rstat,prior=DPM.prior,mcmc=mcmc.conti,state=fit$state,status=FALSE)
    ##calculate cluster
    ss<-fit$state$ss
    num<-fit$state$ncluster
    prop<-table(fit$state$ss)/num
    mu<-fit$state$muclus[1:num]
    sigma<-fit$state$sigmaclus[1:num]
    index<-order(mu,decreasing=T)
    newm<-mu[index]
    news<-sigma[index] 
    newp<-prop[index]
    
    
    if(num==2 | num==1){results[[oo]]<-ss-1}else{
      #hh1<-1/2/sqrt(pi*news[1])+1/2/sqrt(pi*news[2])-2*stats::dnorm(newm[1],newm[2],sd=sqrt(news[1]+news[2]))
      
      #hh2<-1/2/sqrt(pi*news[2])+1/2/sqrt(pi*news[3])-2*stats::dnorm(newm[2],newm[3],sd=sqrt(news[2]+news[3]))
      
      clmatrix<-matrix(rep(0,num*num),ncol=num,nrow=num)
      for(i in 1:num){
        for(j in 1:num){
          clmatrix[i,j]<-stats::dnorm(newm[i],newm[j],sd=sqrt(news[i]+news[j]))
          
        }
      }
      
      
      temp1<-0
      temp2<-0
      com1<-0
      com2<-0
      res1<-0
      res2<-0
      final1<-0
      final2<-0
      latent<-0
      latent[1]<-0
      final1<-c(1,-1,rep(0,num-2))
      final2<-c(0,-1,1,rep(0,num-3))
      
      for(iii in 1:(num-2)){
        dec1<-t(final1)%*%clmatrix%*%final1
        dec2<-t(final2)%*%clmatrix%*%final2
        
        if(dec1>dec2){latent[iii+1]<-1
        }else{
          latent[iii+1]<-0
        }
        if(latent[iii+1]==1){
          ddd<-max(which(latent==0))
          temp1<-rep(0,num)
          temp1[1:ddd]<-1
          com1<-temp1*newp
          com1<-scale(com1,center=F,scale=sum(com1))
          temp1<-rep(0,num)
          temp1[(ddd+1):(iii+1)]<--1
          res1<-temp1*newp
          res1<-scale(res1,center=F,scale=sum(abs(res1)))
          final1<-com1+res1
          temp2<-rep(0,num)
          temp2[(ddd+1):(iii+1)]<--1
          res2<-temp2*newp
          res2<-scale(res2,center=F,scale=sum(abs(res2)))
          com2<-rep(0,num)
          com2[iii+2]<-1
          final2<-com2+res2
        }else{
          temp1<-rep(0,num)
          temp1[1:iii]<-1
          com1<-temp1*newp
          com1<-scale(com1,center=F,scale=sum(com1))
          res1<-rep(0,num)
          res1[iii+1]<--1
          final1<-com1+res1
          final2<-rep(0,num)
          final2[iii+1]<--1
          final2[iii+2]<-1
        }
        
        
      }
      latent[num]<-1
      
      uuu<-max(which(latent==0))
      ss[ss %in% index[(uuu+1):num]]<-0
      ss[ss %in% index[1:uuu]]<-1
      results[[oo]]<-ss
    }
  }
  dpdensitycluster<-do.call(rbind,results)
  mu0=sapply(1:v, function(kk) return(mean(rstat[which(dpdensitycluster[kk,]==0)])))
  mu1=sapply(1:v, function(kk) return(mean(rstat[which(dpdensitycluster[kk,]==1)])))
  var0=sapply(1:v, function(kk) return(stats::var(rstat[which(dpdensitycluster[kk,]==0)])))
  var1=sapply(1:v, function(kk) return(stats::var(rstat[which(dpdensitycluster[kk,]==1)])))
  pro0=sapply(1:v, function(kk) return(mean(rstat[which(dpdensitycluster[kk,]==0)])/length(pvalue)))
  pro1=sapply(1:v, function(kk) return(mean(rstat[which(dpdensitycluster[kk,]==1)])/length(pvalue)))
  for (kk in 1:v){
    if (is.na(var0[kk])){var0[kk]=var1[kk]
    }else if (is.na(var1[kk])){var1[kk]=var0[kk]}
  }
  DPdensityHODC=list()
  DPdensityHODC$mean$mu0=mu0
  DPdensityHODC$mean$mu1=mu1
  DPdensityHODC$variance$var0=var0
  DPdensityHODC$variance$var1=var1
  DPdensityHODC$probility$pro0=pro0
  DPdensityHODC$probility$pro1=pro1
  DPdensityHODC$classification=dpdensitycluster
  return(DPdensityHODC)
}

Try the BANFF package in your browser

Any scripts or data that you put into this service are public.

BANFF documentation built on May 29, 2017, 11:59 a.m.