Nothing
#'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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.