## Not run:
# ## use a few plots to compare densities and hazard functions
#
# ## lognormal
# tSeq <- seq(0.001,0.6,0.001)
# meanlog.true <- -2.4
# sdlog.true <- 0.4
# Yd <- dlnorm(tSeq,meanlog.true,sdlog.true)
# Yh <- hlnorm(tSeq,meanlog.true,sdlog.true)
# max.Yd <- max(Yd)
# max.Yh <- max(Yh)
# Yd <- Yd / max.Yd
# Yh <- Yh / max.Yh
# oldpar <- par(mar=c(5,4,4,4))
# plot(tSeq, Yd, type="n", axes=FALSE, ann=FALSE,
# xlim=c(0,0.6), ylim=c(0,1))
# axis(2,at=seq(0,1,0.2),labels=round(seq(0,1,0.2)*max.Yd,digits=2))
# mtext("Density (1/s)", side=2, line=3)
# axis(1,at=pretty(c(0,0.6)))
# mtext("Time (s)", side=1, line=3)
# axis(4, at=seq(0,1,0.2), labels=round(seq(0,1,0.2)*max.Yh,digits=2))
# mtext("Hazard (1/s)", side=4, line=3, col=2)
# mtext("Lognormal Density and Hazard Functions", side=3, line=2,cex=1.5)
# lines(tSeq,Yd)
# lines(tSeq,Yh,col=2)
# par(oldpar)
#
# ## inverse Gaussian
# tSeq <- seq(0.001,0.6,0.001)
# mu.true <- 0.075
# sigma2.true <- 3
# Yd <- dinvgauss(tSeq,mu.true,sigma2.true)
# Yh <- hinvgauss(tSeq,mu.true,sigma2.true)
# max.Yd <- max(Yd)
# max.Yh <- max(Yh)
# Yd <- Yd / max.Yd
# Yh <- Yh / max.Yh
# oldpar <- par(mar=c(5,4,4,4))
# plot(tSeq, Yd, type="n", axes=FALSE, ann=FALSE,
# xlim=c(0,0.6), ylim=c(0,1))
# axis(2,at=seq(0,1,0.2),labels=round(seq(0,1,0.2)*max.Yd,digits=2))
# mtext("Density (1/s)", side=2, line=3)
# axis(1,at=pretty(c(0,0.6)))
# mtext("Time (s)", side=1, line=3)
# axis(4, at=seq(0,1,0.2), labels=round(seq(0,1,0.2)*max.Yh,digits=2))
# mtext("Hazard (1/s)", side=4, line=3, col=2)
# mtext("Inverse Gaussian Density and Hazard Functions", side=3, line=2,cex=1.5)
# lines(tSeq,Yd)
# lines(tSeq,Yh,col=2)
# par(oldpar)
#
# ## gamma
# tSeq <- seq(0.001,0.6,0.001)
# shape.true <- 6
# scale.true <- 0.012
# Yd <- dgamma(tSeq, shape=shape.true, scale=scale.true)
# Yh <- hgamma(tSeq, shape=shape.true, scale=scale.true)
# max.Yd <- max(Yd)
# max.Yh <- max(Yh)
# Yd <- Yd / max.Yd
# Yh <- Yh / max.Yh
# oldpar <- par(mar=c(5,4,4,4))
# plot(tSeq, Yd, type="n", axes=FALSE, ann=FALSE,
# xlim=c(0,0.6), ylim=c(0,1))
# axis(2,at=seq(0,1,0.2),labels=round(seq(0,1,0.2)*max.Yd,digits=2))
# mtext("Density (1/s)", side=2, line=3)
# axis(1,at=pretty(c(0,0.6)))
# mtext("Time (s)", side=1, line=3)
# axis(4, at=seq(0,1,0.2), labels=round(seq(0,1,0.2)*max.Yh,digits=2))
# mtext("Hazard (1/s)", side=4, line=3, col=2)
# mtext("Gamma Density and Hazard Functions", side=3, line=2,cex=1.5)
# lines(tSeq,Yd)
# lines(tSeq,Yh,col=2)
# par(oldpar)
#
# ## Weibull
# tSeq <- seq(0.001,0.6,0.001)
# shape.true <- 2.5
# scale.true <- 0.085
# Yd <- dweibull(tSeq, shape=shape.true, scale=scale.true)
# Yh <- hweibull(tSeq, shape=shape.true, scale=scale.true)
# max.Yd <- max(Yd)
# max.Yh <- max(Yh)
# Yd <- Yd / max.Yd
# Yh <- Yh / max.Yh
# oldpar <- par(mar=c(5,4,4,4))
# plot(tSeq, Yd, type="n", axes=FALSE, ann=FALSE,
# xlim=c(0,0.6), ylim=c(0,1))
# axis(2,at=seq(0,1,0.2),labels=round(seq(0,1,0.2)*max.Yd,digits=2))
# mtext("Density (1/s)", side=2, line=3)
# axis(1,at=pretty(c(0,0.6)))
# mtext("Time (s)", side=1, line=3)
# axis(4, at=seq(0,1,0.2), labels=round(seq(0,1,0.2)*max.Yh,digits=2))
# mtext("Hazard (1/s)", side=4, line=3, col=2)
# mtext("Weibull Density and Hazard Functions", side=3, line=2,cex=1.5)
# lines(tSeq,Yd)
# lines(tSeq,Yh,col=2)
# par(oldpar)
#
# ## refractory exponential
# tSeq <- seq(0.001,0.6,0.001)
# rate.true <- 20
# rp.true <- 0.01
# Yd <- drexp(tSeq, rate.true, rp.true)
# Yh <- hrexp(tSeq, rate.true, rp.true)
# max.Yd <- max(Yd)
# max.Yh <- max(Yh)
# Yd <- Yd / max.Yd
# Yh <- Yh / max.Yh
# oldpar <- par(mar=c(5,4,4,4))
# plot(tSeq, Yd, type="n", axes=FALSE, ann=FALSE,
# xlim=c(0,0.6), ylim=c(0,1))
# axis(2,at=seq(0,1,0.2),labels=round(seq(0,1,0.2)*max.Yd,digits=2))
# mtext("Density (1/s)", side=2, line=3)
# axis(1,at=pretty(c(0,0.6)))
# mtext("Time (s)", side=1, line=3)
# axis(4, at=seq(0,1,0.2), labels=round(seq(0,1,0.2)*max.Yh,digits=2))
# mtext("Hazard (1/s)", side=4, line=3, col=2)
# mtext("Refractory Exponential Density and Hazard Functions", side=3, line=2,cex=1.5)
# lines(tSeq,Yd)
# lines(tSeq,Yh,col=2)
# par(oldpar)
#
# ## log logistic
# tSeq <- seq(0.001,0.6,0.001)
# location.true <- -2.7
# scale.true <- 0.025
# Yd <- dllogis(tSeq, location.true, scale.true)
# Yh <- hllogis(tSeq, location.true, scale.true)
# max.Yd <- max(Yd)
# max.Yh <- max(Yh)
# Yd <- Yd / max.Yd
# Yh <- Yh / max.Yh
# oldpar <- par(mar=c(5,4,4,4))
# plot(tSeq, Yd, type="n", axes=FALSE, ann=FALSE,
# xlim=c(0,0.6), ylim=c(0,1))
# axis(2,at=seq(0,1,0.2),labels=round(seq(0,1,0.2)*max.Yd,digits=2))
# mtext("Density (1/s)", side=2, line=3)
# axis(1,at=pretty(c(0,0.6)))
# mtext("Time (s)", side=1, line=3)
# axis(4, at=seq(0,1,0.2), labels=round(seq(0,1,0.2)*max.Yh,digits=2))
# mtext("Hazard (1/s)", side=4, line=3, col=2)
# mtext("Log Logistic Density and Hazard Functions", side=3, line=2,cex=1.5)
# lines(tSeq,Yd)
# lines(tSeq,Yh,col=2)
# par(oldpar)
# ## End(Not run)
Run the code above in your browser using DataLab