# Simulate Weibull distribution 

# Parameters:
# 1. shape -- alpha  = 2
# 2. scale -- lambda = 1/6 (in R notations scale =6)


y <- seq(0,20,by = 0.001)

par(mfcol=c(1,3))
plot(y,dweibull(y,2,6),type="l",main = "Weibull density",ylab="density",xlab="hours")
plot(y,pweibull(y,2,6),type="l",main = "Weibull cdf",ylab="Probability",xlab="hours")
plot(y,1-pweibull(y,2,6),type="l",main = "Weibull survival",ylab="Probability",xlab="hours")


#  1. Probability of 5 to 6 hour survival

pweibull(6,2,6) - pweibull(5,2,6)

par(mfcol=c(1,1))
plot(y,1-pweibull(y,2,6),type="l",main = "Weibull survival",ylab="Probability",xlab="hours",lwd=2)
segments(5,0,5,1-pweibull(5,2,6),lty=2)
segments(6,0,6,1-pweibull(6,2,6),lty=2)

segments(5,1-pweibull(6,2,6),6,1-pweibull(6,2,6),lty=2)
arrows(5,1-pweibull(6,2,6),5,1-pweibull(5,2,6),col=2,length=.1,code=3,lwd=2)

plot(y,1-pweibull(y,2,6),type="l",main = "Weibull survival",ylab="Probability",xlab="hours",lwd=2,xlim=c(0,10),ylim=c(.2,.6))
segments(5,0,5,1-pweibull(5,2,6),lty=2)
segments(6,0,6,1-pweibull(6,2,6),lty=2)
text(5.5,.5,paste(round(1-pweibull(5,2,6),3)),col=3)
text(6.5,1-pweibull(6,2,6),paste(round(1-pweibull(6,2,6),3)),col=3)
segments(5,1-pweibull(6,2,6),6,1-pweibull(6,2,6),lty=2)
arrows(5,1-pweibull(6,2,6),5,1-pweibull(5,2,6),col=2,length=.1,code=3,lwd=2)

plot(y,dweibull(y,2,6),type="l",main = "Weibull density",ylab="",xlab="hours",lwd=2)
x5.0 <- c(5,seq(5,6,by=0.01),6,5)
y5.0 <- c(0,dweibull(seq(5,6,by=0.01),2,6),0,0)
polygon(x5.0,y5.0,density=20,col=3)
text(7.5,.13,paste(round(pweibull(6,2,6)-pweibull(5,2,6),4)),col=3,cex=1.2)

x5.1 <- c(5,5,6,6,5)
y5.1 <- c(0,dweibull(c(5.5,5.5),2,6),0,0)
polygon(x5.1,y5.1,angle=135,density=20,col=2)
text(7.5,.125,paste(round(dweibull(5.5,2,6),4)),col=2,cex=1.2)


#  2. Probability of 300 to 360 minute survival

plot(y*60,1-pweibull(y,2,6),type="l",main = "Weibull survival",ylab="Probability",xlab="minutes",lwd=2)
segments(300,0,300,1-pweibull(5,2,6),lty=2)
segments(360,0,360,1-pweibull(6,2,6),lty=2)
segments(300,1-pweibull(6,2,6),360,1-pweibull(6,2,6),lty=2)
arrows(300,1-pweibull(6,2,6),300,1-pweibull(5,2,6),col=2,length=.1,code=3,lwd=2)

plot(y*60,dweibull(y,2,6)/60,type="l",main = "Weibull density",ylab="",xlab="hours",lwd=2)
x5.1 <- c(300,300,360,360,300)
y5.1 <- c(0,dweibull(c(5.5,5.5),2,6)/60,0,0)
polygon(x5.1,y5.1,density=40,col=2)
text(700,.0021,paste(round(dweibull(5.5,2,6),4)," = 60 *",round(dweibull(5.5,2,6)/60,4)),col=2,cex=1.2)

#  3. hazard 


# Weibull hazard: 
# h(t) = lambda * alpha * ( lambda * t) ^ (alpha - 1)  
#      = 1/6 * 2 * 1/6 * t =  1/18 * t

plot(y,dweibull(y,2,6) / (1-pweibull(y,2,6)),type="l",main = "Weibull hazard",ylab="",xlab="hours",lwd=2)
segments(5,0,5,dweibull(5,2,6) / (1-pweibull(5,2,6)),lty=2)
text(8,.25,paste(round(dweibull(5,2,6) / (1-pweibull(5,2,6)),3)," = 5 / 18" ),col=3,cex=1.2)

plot(y,dweibull(y,2,6),type="l",ylab="",xlab="hours",lwd=2,main = "Computation of h(5) = 0.278 ")

x5.0 <- c(5,seq(5,20,by=0.01),20,5)
y5.0 <- c(0,dweibull(seq(5,20,by=0.01),2,6),0,0)
polygon(x5.0,y5.0,density=20,col=3)

x5.1 <- c(5,seq(5,6,by=0.01),6,5)
y5.1 <- c(0,dweibull(seq(5,6,by=0.01),2,6),0,0)
polygon(x5.1,y5.1,angle=135,density=20,col=2)
text(7.5,.13,paste(round( (pweibull(6,2,6)-pweibull(5,2,6)) / ( 1 - pweibull(5,2,6)),4)),col=3,cex=1.2)


plot(y,dweibull(y,2,6),type="l",ylab="",xlab="hours",lwd=2,main = "Better computation of h(5) = 0.2778 ")

x5.0 <- c(5,seq(5,20,by=0.01),20,5)
y5.0 <- c(0,dweibull(seq(5,20,by=0.01),2,6),0,0)
polygon(x5.0,y5.0,density=20,col=3)

x5.1 <- c(5,seq(5,5.1,by=0.01),5.1,5)
y5.1 <- c(0,dweibull(seq(5,5.1,by=0.01),2,6),0,0)
polygon(x5.1,y5.1,angle=135,density=20,col=2)
text(7.5,.13,paste(round( (pweibull(5.1,2,6)-pweibull(5,2,6)) / ( 1 - pweibull(5,2,6)),5)),col=2,cex=1.2)

segments(6,.129,9,.129,col=2) ; text(7.5,.138,paste(round( 10*(pweibull(5.1,2,6)-pweibull(5,2,6)) / ( 1 - pweibull(5,2,6)),4)),col=3,cex=1.2)


# Weibull hazard in minutes

plot(y*60,dweibull(y,2,6)/60,type="l",main = paste("Weibull hazard in minutes  - h(300) = ",
round( (dweibull(5,2,6) / 60) / ( 1 - pweibull(5,2,6)),7)),ylab="density",xlab="minutes",lwd=2)

x5.0 <- c(5,seq(5,20,by=0.01),20,5)*60
y5.0 <- c(0,dweibull(seq(5,20,by=0.01),2,6),0,0) / 60
polygon(x5.0,y5.0,density=20,col=3)

x5.1 <- c(300,seq(300,301,by=0.01),301,300)
y5.1 <- c(0,dweibull(seq(300,301,by=0.01)/60,2,6),0,0) / 60
polygon(x5.1,y5.1,angle=135,density=20,col=2)

text(500,0.002,paste(round( (pweibull(5+1/60,2,6)-pweibull(5,2,6))  / ( 1 - pweibull(5,2,6)),7)),col=2,cex=1.2)