rm(list=ls())
setwd("/Users/liliwang/Box/Cusum_simulation/restart3/For_paper/results/head_starts")
#date="20180311"
rho_list=seq(0,0.95,0.05)
mu_list=c(log(3),log(2),log(1.5),log(1.2),log(1),log(1),-log(1.2),-log(1.5),-log(2),-log(3)) #1-10
mu_list_str=c(paste0("log",c(3,2,1.5,1.2,1)),paste0("-log",c(1,1.2,1.5,2,3)))
data=NULL
for(jobid in 1:10){
data<-rbind(data,cbind(read.csv(paste0("rho_res",jobid,".csv"))[,2:4],rep(mu_list_str[jobid],length(rho_list))))
}
colnames(data)<-c("rho","ARL","alpha","mu")
data$ARL=data$ARL/365
library(ggplot2)
data_up=data[1:100,]
data_up$mu<-factor(data_up$mu,levels=paste0("log",c(3,2,1.5,1.2,1)))
Sim2_up_ARL<-ggplot(data_up,aes(x=rho,y=ARL,shape=mu))+geom_line(size=0.3)+geom_point()+
ggtitle("theta: log2")+ylim(range(data$ARL))+
theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))
Sim2_up_ARL
ggsave(filename="Sim2_11.pdf", plot=Sim2_up_ARL,dpi = 300, unit="in", width=6,height=6)
Sim2_up_alpha<-ggplot(data_up,aes(x=rho,y=alpha,shape=mu))+geom_line(size=0.3)+geom_point()+
ggtitle("theta: log2")+ylim(range(data$alpha))+
theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))
Sim2_up_alpha
ggsave(filename="Sim2_21.pdf", plot=Sim2_up_alpha,dpi = 300, unit="in", width=6,height=6)
data_low=data[101:200,]
data_low$mu<-factor(data_low$mu,levels=paste0("-log",c(3,2,1.5,1.2,1),""))
Sim2_low_ARL<-ggplot(data_low,aes(x=rho,y=ARL,shape=mu))+geom_line(size=0.3)+geom_point()+
ggtitle("theta: -log2")+ylim(range(data$ARL))+
theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))
Sim2_low_ARL
ggsave(filename="Sim2_12.pdf", plot=Sim2_low_ARL,dpi = 300, unit="in", width=6,height=6)
Sim2_low_alpha<-ggplot(data_low,aes(x=rho,y=alpha,shape=mu))+geom_line(size=0.3)+geom_point()+
ggtitle("theta: -log2")+ylim(range(data$alpha))+
theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))
Sim2_low_alpha
ggsave(filename="Sim2_22.pdf", plot=Sim2_low_alpha,dpi = 300, unit="in", width=6,height=6)
rm(list=ls()) #remove list
setwd("/Users/liliwang/Box/Cusum_simulation/restart3/For_paper/results/control_limits")
theta=log(c(3,2,1.5,1.2))
theta=cbind(theta,-theta)
theta_list=c("log3","log2","log1.5","log1.2")
Nevent=seq(20,200,10)*0.1
mu=log(1)
theta_star<-data<-NULL
for(jobid in 1:4){
data<-rbind(data,cbind(Nevent,data.frame(read.csv(paste0("selected_h_test",jobid,"_2.csv"),header=T)[,3:4]),rep(theta_list[jobid],length(Nevent)),rep(theta[jobid],length(Nevent))))
#theta_star<-rbind(theta_star,cbind(
#rep(log((exp(theta[jobid,1])-exp(mu))/(theta[jobid,1]-mu)),length(Nevent)),
#rep(log((exp(theta[jobid,2])-exp(mu))/(theta[jobid,2]-mu)),length(Nevent))))
}
library(ggplot2)
rownames(data)<-NULL
colnames(data)[1:5]<-c("N_event","h_up","h_low","theta","adj")
data=data.frame(data)
### Change color of both line and points
# Change line type and point type, and use thicker line and larger points
# Change points to circles with white fill
#ggplot(data=dat, aes(x=time, y=total_bill, group=1)) +
# geom_line(colour="red", linetype="dashed", size=1.5) +
# geom_point(colour="red", size=4, shape=21, fill="white")
pl_up<-ggplot(data,aes(x=N_event,y=h_up,shape=theta))+geom_line(size=0.7)+xlab("Expected 1-year failtures")+geom_point(size=2,fill="white")+
ylab("h values") + ylim(c(range(c(data$h_up,data$h_low)))) +
ggtitle("Worse than expected")+theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))+
scale_shape_manual(name="theta",values=21:24,label=theta_list)
#theme(text = element_text(size=20),
# axis.text.x = element_text(angle=90, hjust=1))
ggsave(filename="Fig1_1.pdf", plot=pl_up,dpi = 300, unit="in", width=6,height=6 )
#data$theta=paste0("-",data$theta)
pl_low<-ggplot(data,aes(x=N_event,y=h_low,shape=theta))+geom_line(size=0.7)+
geom_point(size=2,fill="white")+xlab("Expected 1-year failtures") +
ylab("h values") + ylim(c(range(c(data$h_up,data$h_low)))) +
ggtitle("Better than expected")+theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))+
scale_shape_manual(name="theta",values=21:24,label=paste0("-",theta_list))
# axis.text.x = element_text(angle=90, hjust=1))
ggsave(filename="Fig1_2.pdf", plot=pl_low,dpi = 300, unit="in", width=6,height=6)
data$up_adj=data$h_up*data$adj
data$low_adj=data$h_low*data$adj
pl_up2<-ggplot(data,aes(x=N_event,y=up_adj,shape=theta))+geom_line(size=0.7)+xlab("Expected 1-year failtures")+geom_point(size=2,fill="white")+
ylab("L values") + ylim(c(range(c(data$up_adj,data$low_adj)))) +
ggtitle("Worse than expected")+theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))+
scale_shape_manual(name="theta",values=21:24,label=theta_list)
#theme(text = element_text(size=20),
# axis.text.x = element_text(angle=90, hjust=1))
ggsave(filename="Fig2_1.pdf", plot=pl_up2,dpi = 300, unit="in", width=6,height=6 )
#data$theta=paste0("-",data$theta)
pl_low2<-ggplot(data,aes(x=N_event,y=low_adj,shape=theta))+geom_line(size=0.7)+
geom_point(size=2,fill="white")+xlab("Expected 1-year failtures") +
ylab("L values") + ylim(c(range(c(data$up_adj,data$low_adj)))) +
ggtitle("Better than expected")+theme(plot.title = element_text(hjust = 0.5,size=20),
legend.title=element_text(size=18),
legend.text=element_text(size=15))+
scale_shape_manual(name="theta",values=21:24,label=paste0("-",theta_list))
# axis.text.x = element_text(angle=90, hjust=1))
ggsave(filename="Fig2_2.pdf", plot=pl_low2,dpi = 300, unit="in", width=6,height=6)
##Check whether the h values have satisfy: h/(exp(mu)-exo(theta_start)) to be similar
#data
#theta_star
#dim(theta_star)
#dim(data)
#(data$flag_up_h/abs(1-exp(theta_star[,1])))
#(data$flag_low_h/abs(1-exp(theta_star[,2])))
