demo/sim2026_runCode.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

# Simulation Run File
#     Change Run number to change sample sizes and proportion censored
#setwd("c:/R/work/simBPCP")
#getwd()

#Run_number<- 2
#NSIM<- 3
#PARMTYPE<-"difference"
#PARMTYPE<-"efflogs"

if (PARMTYPE=="difference"){
  BETAEQUAL<- 0
  BETAFunc<-function(S1,S2){ S2-S1 }
} else if (PARMTYPE=="efflogs"){
  BETAEQUAL<- 0
  BETAFunc<-function(S1,S2){ 
    beta<-1-log(S2)/log(S1) 
    beta[S1==1 & S2==1]<- 0
    beta
  }
} else stop("rewrite program, give BETAFUNC for other PARMTYPES")


# First simulation that took more than an overnight to run 
# had two errors (D15 and E18) had the error that two random 
# times where about 1e-11 apart that gives an error from the 
# software. I changed that error to a warning and reran 
# those two simulation runs (D15 and E18). 


source("Sim2026_functions.R")
#library(NoSleepR)
#nosleep_on()



DO_MIDP<-TRUE


Run_parameters<-matrix(c(
  30,20,40,30,20,40,90,60,120,90,60,120,150,100,200,150,100,200,
  30,40,20,30,40,20,90,120,60,90,120,60,150,200,100,150,200,100,
  rep(.40,3),rep(.10,3),rep(.40,3),rep(.10,3),rep(.40,3),rep(.10,3)),18,3,byrow=FALSE,
  dimnames=list(paste0("Run=", 1:18),c("N1","N2","pCENS")))

Run_i<-Run_parameters[Run_number,]

N1<- Run_i["N1"]
N2<- Run_i["N2"]
pCENS<- Run_i["pCENS"]


##############################
## Simulation for Scenario 1
##############################
t0<-proc.time()

TTIME<- 0.8

outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen1,funcs=f1,domidp=DO_MIDP)
betaTRUE<- BETAFunc(f1$S1(TTIME), f1$S2(TTIME) )
s1n30t08<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.0
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen1,funcs=f1,domidp=DO_MIDP)
betaTRUE<- BETAFunc(f1$S1(TTIME), f1$S2(TTIME) )
s1n30t10<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.2
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen1,funcs=f1,domidp=DO_MIDP)
betaTRUE<- BETAFunc(f1$S1(TTIME), f1$S2(TTIME) )
s1n30t12<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.4
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen1,funcs=f1,domidp=DO_MIDP)
betaTRUE<- BETAFunc(f1$S1(TTIME), f1$S2(TTIME) )
s1n30t14<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.8
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen1,funcs=f1,domidp=DO_MIDP)
betaTRUE<- BETAFunc(f1$S1(TTIME), f1$S2(TTIME) )
s1n30t18<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))

results_Sc1<-rbind(s1n30t08,
      s1n30t10,
      s1n30t12,
      s1n30t14,
      s1n30t18)


t1<-proc.time()
t1-t0



##############################
## Simulation for Scenario 2
##############################


set.seed(312*Run_number)

TTIME<- 0.8

outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen2,funcs=f2,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f2$S1(TTIME),f2$S2(TTIME))
s2n30t08<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.0
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen2,funcs=f2,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f2$S1(TTIME),f2$S2(TTIME))
s2n30t10<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.2
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen2,funcs=f2,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f2$S1(TTIME),f2$S2(TTIME))
s2n30t12<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.4
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen2,funcs=f2,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f2$S1(TTIME),f2$S2(TTIME))
s2n30t14<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.8
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen2,funcs=f2,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f2$S1(TTIME),f2$S2(TTIME))

s2n30t18<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))

results_Sc2<-rbind(s2n30t08,
      s2n30t10,
      s2n30t12,
      s2n30t14,
      s2n30t18)

t2<-proc.time()
t2-t0



##############################
## Simulation for Scenario 3
##############################

set.seed(123*Run_number)
TTIME<- 0.8

outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen3,funcs=f3,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f3$S1(TTIME),f3$S2(TTIME))
s3n30t08<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.0
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen3,funcs=f3,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f3$S1(TTIME),f3$S2(TTIME))
s3n30t10<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.2
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen3,funcs=f3,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f3$S1(TTIME),f3$S2(TTIME))
s3n30t12<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.4
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen3,funcs=f3,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f3$S1(TTIME),f3$S2(TTIME))
s3n30t14<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.8
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen3,funcs=f3,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f3$S1(TTIME),f3$S2(TTIME))
s3n30t18<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))

results_Sc3<-rbind(s3n30t08,
      s3n30t10,
      s3n30t12,
      s3n30t14,
      s3n30t18)
t3<-proc.time()
t3-t0


##############################
## Simulation for Scenario 4
##############################


set.seed(912*Run_number)
TTIME<- 0.8

outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen4,funcs=f4,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f4$S1(TTIME),f4$S2(TTIME))
s4n30t08<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.0
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen4,funcs=f4,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f4$S1(TTIME),f4$S2(TTIME))
s4n30t10<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.2
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen4,funcs=f4,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f4$S1(TTIME),f4$S2(TTIME))
s4n30t12<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.4
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen4,funcs=f4,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f4$S1(TTIME),f4$S2(TTIME))
s4n30t14<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))


TTIME<- 1.8
outTable<-sim(nsim=NSIM,TESTTIME=TTIME,n1=N1,n2=N2,
              pCens=pCENS,simData=simScen4,funcs=f4,domidp=DO_MIDP)
betaTRUE<-BETAFunc(f4$S1(TTIME),f4$S2(TTIME))
s4n30t18<-c(t=TTIME,beta=betaTRUE,
            simSummary(outTable,beta=betaTRUE))

results_Sc4<-rbind(s4n30t08,
      s4n30t10,
      s4n30t12,
      s4n30t14,
      s4n30t18)
t4<-proc.time()


Results<- rbind(results_Sc1,
                results_Sc2,
                results_Sc3,
                results_Sc4)
Results


save(Results,file=paste0("simResults/Results_",PARMTYPE,"_Run_",Run_number,".RData"))
#load(file=paste0("simResults/Results_",PARMTYPE,"_Run_",Run_number,".RData"))

t4-t0

#nosleep_off()