# 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