# for bic simulations
x.quad.0 <- matrix(scan("K:/srcis2/forward.r/simulation/x.quad.0.txt",0),ncol=252,byrow=TRUE)
x.quad.70 <- matrix(scan("K:/srcis2/forward.r/simulation/x.quad.70.txt",0),ncol=252,byrow=TRUE)
xmatrix.0<-x.quad.0[,1:21]                                # 150 by 21
xmatrix.70<-x.quad.70[,1:21]                              # rho=.70
hbeta<-matrix(0,5,14)                                     # 5 rows of 14 cols
hbeta[1,]<-c(rep(-1,14))
hbeta[2,]<-c(7,14,rep(-1,12))                             # true betas with
hbeta[3,]<-c(6,7,8,13,14,15,rep(-1,8))                    # -1's to make same length
hbeta[4,]<-c(5,6,7,8,9,12,13,14,15,16,rep(-1,4))
hbeta[5,]<-c(4,5,6,7,8,9,10,11,12,13,14,15,16,17)
res<-matrix(0,5,14)
me.out<-matrix(rep(0,400),nrow=100,ncol=4,byrow=T)
me.out.21.0<-array(0,dim=c(100,4,5))
me.out.21.70<-array(0,dim=c(100,4,5))
jac<-rep(0,100)
drv <-"K:\\srcis2\\forward.r\\simulation\\h0_0.rs35.txt"
 h <- c("0","1","2","3","4")
for(m in 1:5){                                            # loop through data sets
substr(drv,33,33)<-h[m]                                   # puts in 1,2,3, or 4
read.table(drv,header=T)->hdata                           # rep,mu,y
matrix(hdata[,3],nrow=100,ncol=150,byrow=T)->hdata.y      # each row is a y data set
mu<-hdata[1:150,2]                                        # true means
beta<-hbeta[m,]                                           # indicator of nonzero beta

for (i in 1:100){                                         # loop through data sets
out.bic<-bic.sim(x=xmatrix.0,y=hdata.y[i,])
st<-out.bic$size                                          # number of fitted x's
cs<-length(intersect(out.bic$x.ind,beta))                 # number correct
fs<-st-cs                                                 # number false
me.out[i,1]<-mean((out.bic$mod$fitted-mu)^2)              # model error
me.out[i,2]<-st
me.out[i,3]<-fs
me.out[i,4]<-sum(out.bic$mod$resid^2)/out.bic$mod$df.residual
}       # ends i loop
me.out.21.0[,,m]<-round(me.out,4)

mod.size<-sum(beta != -1)                             # true model size
for (i in 1:100){
if(me.out[i,3]+mod.size > 0) jac[i]<-(me.out[i,2]-me.out[i,3])/(me.out[i,3]+mod.size) else jac[i]<-1}
res[m,1]<-0                                           # rho
res[m,2]<-m-1                                         # model h0, h1, ..
res[m,3]<-mean(me.out[,2])                            # model size
res[m,4]<-sd(me.out[,2])/10                           # se of the mean
res[m,5]<-mean(me.out[,1])                            # me
res[m,6]<-sd(me.out[,1])/10
res[m,7]<-mean(me.out[,3]/(1+me.out[,2]))             # fsr as mean of fsr's
res[m,8]<-sd(me.out[,3]/(1+me.out[,2]))/10
res[m,9]<-mean(me.out[,4])                            # mse of selected model
res[m,10]<-sd(me.out[,4])/10
if(mod.size > 0) csr<-(me.out[,2]-me.out[,3])/mod.size else csr<-rep(1,100)
res[m,11]<-mean(csr)
res[m,12]<-sd(csr)/10
res[m,13]<-mean(jac)
res[m,14]<-sd(jac)/10
}                                                      # for m loop of 5 models
bic.cal.35.0<-data.frame(rho=res[,1],H=res[,2],size=res[,3],
me=res[,5],fsr_mr=res[,7],mse=res[,9],csr=res[,11],jac=res[,13],
se.size=res[,4],se.me=res[,6],se.fsr=res[,8],se.mse=res[,10],se.csr=res[,12],se.jac=res[,14])
print(round(bic.cal.35.0,3))

# Now for rho=.70  Note changes in drv, and xmatrix.0 to xmatrix.70
drv <-"K:\\srcis2\\forward.r\\simulation\\h0_70.rs35.txt"
for(m in 1:5){                                            # loop through data sets
substr(drv,33,33)<-h[m]                                   # puts in 1,2,3, or 4
read.table(drv,header=T)->hdata                           # rep,mu,y
matrix(hdata[,3],nrow=100,ncol=150,byrow=T)->hdata.y      # each row is a y data set
mu<-hdata[1:150,2]                                        # true means
beta<-hbeta[m,]                                           # indicator of nonzero beta

for (i in 1:100){
out.bic<-bic.sim(x=xmatrix.70,y=hdata.y[i,])
st<-out.bic$size                                          # number of fitted x's
cs<-length(intersect(out.bic$x.ind,beta))                 # number correct
fs<-st-cs                                                 # number false
me.out[i,1]<-mean((out.bic$mod$fitted-mu)^2)              # model error
me.out[i,2]<-st
me.out[i,3]<-fs
me.out[i,4]<-sum(out.bic$mod$resid^2)/out.bic$mod$df.residual
}       # ends i loop
me.out.21.70[,,m]<-round(me.out,4)

mod.size<-sum(beta != -1)                             # true model size
for (i in 1:100){
if(me.out[i,3]+mod.size > 0) jac[i]<-(me.out[i,2]-me.out[i,3])/(me.out[i,3]+mod.size) else jac[i]<-1}
res[m,1]<-0.7                                           # rho
res[m,2]<-m-1                                         # model h0, h1, ..
res[m,3]<-mean(me.out[,2])                            # model size
res[m,4]<-sd(me.out[,2])/10                           # se of the mean
res[m,5]<-mean(me.out[,1])                            # me
res[m,6]<-sd(me.out[,1])/10
res[m,7]<-mean(me.out[,3]/(1+me.out[,2]))             # fsr as mean of fsr's
res[m,8]<-sd(me.out[,3]/(1+me.out[,2]))/10
res[m,9]<-mean(me.out[,4])                            # mse of selected model
res[m,10]<-sd(me.out[,4])/10
if(mod.size > 0) csr<-(me.out[,2]-me.out[,3])/mod.size else csr<-rep(1,100)
res[m,11]<-mean(csr)
res[m,12]<-sd(csr)/10
res[m,13]<-mean(jac)
res[m,14]<-sd(jac)/10
}                                                      # for m loop of 5 models
bic.cal.35.70<-data.frame(rho=res[,1],H=res[,2],size=res[,3],
me=res[,5],fsr_mr=res[,7],mse=res[,9],csr=res[,11],jac=res[,13],
se.size=res[,4],se.me=res[,6],se.fsr=res[,8],se.mse=res[,10],se.csr=res[,12],se.jac=res[,14])
print(round(bic.cal.35.70,3))