# This document contains four main files: # function 'rdata.hmm' (rdata1.hmm) is the HMM data generator # function 'bwfw.hmm' (bwfw1.hmm) realizes the forward-backward procedure # function 'em.hmm' (em1.hmm) realizes the E-M algorithm # function 'mt.hmm' realizes the new multiple testing procedures # (given the estimates of HMM parameters) # Examples on how to use these functions are given below. # R codes for other multiple testing procedures are also included for the comprehensive Example 5: # function 'mt.bh' realizes the BH step-up procedure # function 'mt.gw' realizes the adaptive p-value procedure # More detailed instructions are given in each separate file. ############################### ###### EXAMPLES ########### ############################### source("rdata1.hmm.R.txt") source("rdata.hmm.R.txt") source("bwfw1.hmm.R.txt") source("bwfw.hmm.R.txt") source("em1.hmm.R.txt") source("em.hmm.R.txt") source("mt.hmm.R.txt") source("mt.gw.R.txt") source("mt.bh.R.txt") # the number of observations NUM1<-1000 NUM<-3000 # the prespecified FDR level q=0.10 #################################### # Example 1: the HMM data generator #################################### ## (1.a) one-component alternative # the initial state distribution pii<-c(1, 0) # the transition matrix A<-matrix(c(0.95, 0.05, 0.2, 0.80), 2, 2, byrow=T) # the null distribution f0<-c(0, 1) # the alternative distribution f1<-c(1, 1) # the HMM data set.seed(123456) rdata1<-rdata.hmm(NUM1, pii, A, f0, 1, f1) # the observed values x1<-rdata1\$o # the unobserved states theta1<-rdata1\$s ## (1.b) three-component alternative # the initial state distribution pii<-c(1, 0) # the transition matrix A<-matrix(c(0.95, 0.05, 0.2, 0.80), 2, 2, byrow=T) # the null distribution f0<-c(0, 1) # the alternative distribution pc<-c(0.4, 0.3, 0.3) f1.v<-matrix(c(-1, 1, 2, 1, 3, 1), 3, 2, byrow=TRUE) # the HMM data rdata2<-rdata.hmm(NUM, pii, A, f0, pc, f1.v) ############################################ # Example 2: the forward-backward procedure ############################################ ## (2.a) one-component alternative x1<-rdata1\$o fb.res1<-bwfw1.hmm(x1, pii, A, f0, f1) # the backward variable backward.var<-fb.res1\$bw # the backward variable forward.var<-fb.res1\$fw # the LSI variable lsi.var<-fb.res1\$lsi ## (2.b) three-component alternative x2<-rdata2\$o fb.res2<-bwfw.hmm(x2, pii, A, pc, f0, f1.v) ################################################### # Example 3: the E-M Algorithm in a normal mixture ################################################### ## (3.a) one-component alternative L<-1 # the EM algorithm em.res1<-em1.hmm(x1, maxiter=500) # the estimates for HMM parameters em.res1\$A em.res1\$f1 # the number of interations em.res1\$n1 ## (3.b) three-component alternative L<-3 em.res2<-em.hmm(x2, L=3, maxiter=500) ############################################ # Example 4: The AP, OR and PI procedures ############################################ # the p-values pv1<-1-pnorm(x1, 0, 1) # the LSI values ## (4.a) the AP procedure A0<-em.res1\$A p0<-A0[2, 1]/(A0[1, 2]+A0[2, 1]) gw.res<-mt.gw(pv1, q, p0) # the decision rule, gw.de<-gw.res\$de ## (4.b) the OR procedure lsi.or<-fb.res1\$lsi or.res<-mt.hmm(lsi.or, q) # the decision rule or.de<-or.res\$de ## (4.c) the PI procedure lsi.pi<-em.res1\$lf pi.res<-mt.hmm(lsi.pi, q) # the decision rule pi.de<-pi.res\$de ########################################## # Example 5: the analysis of the HIV data ########################################## # the FDR level q<-0.05 # the HIV data hiv<-scan("hivdata.txt") # the p-values hiv.NUM<-length(hiv) hiv.pv<-1-pnorm(hiv, 0, 1) for (j in 1:hiv.NUM) { hiv.pv[j]<-2*min(hiv.pv[j], (1-hiv.pv[j])) } # the LSI values hiv.res<-em.hmm(hiv, L=2, maxiter=500) hiv.lsi<-hiv.res\$lf # estimates of the proportion of non-nulls A0.hiv<-hiv.res\$A p0.hiv<-A0.hiv[2, 1]/(A0.hiv[1, 2]+A0.hiv[2, 1]) # the BH procedure hiv.bh<-mt.bh(hiv.pv, q) # the threshold for p-value hiv.bh\$th # the number of rejections hiv.bh\$nr # the indices of rejected hypotheses hiv.bh\$re # the AP procedure hiv.ap<-mt.gw(hiv.pv, q, p0.hiv) # the threshold for p-value hiv.pi\$th # the number of rejections hiv.ap\$nr # the indices of rejected hypotheses hiv.ap\$re # the PI procedure hiv.pi<-mt.hmm(hiv.lsi, q) # the threshold for lsi hiv.pi\$th # the number of rejections hiv.pi\$nr # the indices of rejected hypotheses hiv.pi\$re