### 4.3.1 ###

 source("ptest2s.R")
 set.seed(0)
 B <- 5000
 T <- ptest2s(Y,X,B,"Student")

 maxT <- apply(T,1,max)
 sum(maxT[-1]>=maxT[1])/(B-1)	#p-value sumT

 sumT <- apply(T,1,sum)
 sum(sumT[-1]>=sumT[1])/(B-1)	#p-value maxT

 ssT2 <- apply(T,1,function(row){sum(sign(row)*(row^2))})
 sum(ssT2[-1]>= ssT2[1])/(B-1)	#p-value ssT2


 alpha <- 0.05
 c1 <- sort(T[-1,1])[B*(1-alpha)]
 c2 <- sort(T[-1,2])[B*(1-alpha)]
 T[1,] > c(c1,c2)


 ind <- (T[,1]>c1|T[,2]>c2)
 sum(sumT[ind]>=sumT[1])/(B-1)	#p-value csumT

 sum(ssT2[ind]>=ssT2[1])/(B-1)	#p-value cssT2


 eps <- 0.5
 ct <- qt(1-alpha,length(Y)-2)
 delta <- eps/(mean(apply(Y,2,sd)))
 T[1,] + delta > c(ct,ct)

 ind2 <- (T[,1]+delta>ct & T[,2]+delta>ct)
 sum(maxT[ind2]>=maxT[1])/(B-1)	# p-value

### 4.4.1 ###

 load("rats.Rdata")
 source("pstocRs.R") 
 set.seed(0)
 B <- 5000
 T_hi <- pstocRs(X,Y,B)

 load("t2p.R")
 P_hi <- T_hi
 r = length(unique(X))
 for (i in 1:(r-1)){ P_hi[,,i] <- t2p(T_hi[,,i]) }

 T_h <- P_hi[,,1]
 T_h <- apply(P_hi,c(1,2),function(x){-2*log(prod(x))})
 P_h<-t2p(T_h)
 P_h[1,]

 T <- P_hi[,1,1]
 T <- apply(P_h,1,function(x){-2*log(prod(x))})
 sum(T[-1]>=T[1])/(B-1)