############################################### ## Author: Joshua M. Tebbs ## Date: 11 November 2018 ## Update: 30 November 2018 ## STAT 512 course notes: R Code ## Chapter 6 ############################################### # Figure 6.1 # Page 3 # Left y = seq(0,1,0.01) pdf = rep(1,length(y)) plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxt='n',yaxt='n',xlim=c(0,1.5),ylim=c(0,1.5),cex.lab=1.25) axis(side=1,at=c(0,1),labels=c("0","1")) axis(side=2,at=c(0,1),labels=c("0","1")) lines(c(1,1.5),c(0,0),lty=1) # Right y = seq(0.001,1,0.001) h = -log(y) plot(y,h,type="l",xlab="y",ylab="h(y)",xlim=c(0,1),ylim=c(0,6),xaxt='n',cex.lab=1.25) axis(side=1,at=c(0,0.5,1),labels=c("0","0.5","1")) abline(h=0) # Figure 6.2 # Page 5 # Left y = seq(0,6,0.01) pdf = dexp(y,1) plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxt='n',yaxt='n',xlim=c(0,6),ylim=c(0,0.9),cex.lab=1.25) axis(side=1,at=c(0),labels=c("0")) axis(side=2,at=c(0),labels=c("0")) abline(h=0) abline(v=0,lty=2) # Right y = seq(0,4,0.001) m = 3 h = y^(1/m) plot(y,h,type="l",xlab="y",ylab="h(y)",xlim=c(0,4),ylim=c(0,1.8),xaxt='n',yaxt='n',cex.lab=1.25) axis(side=1,at=c(0),labels=c("0")) axis(side=2,at=c(0),labels=c("0")) abline(h=0) # Figure 6.3 # Page 6 y = seq(0,20,0.01) plot(y,dweibull(y,shape=1.5,scale=5^(1/1.5)),type="l",lty=1,xlab="u",ylab="PDF",cex.lab=1.25) lines(y,dweibull(y,shape=2,scale=50^(1/2)),lty=4) lines(y,dweibull(y,shape=2.5,scale=100^(1/2.5)),lty=8) abline(h=0) legend(10,0.20,lty=c(1,4,8),c(expression(paste(m,"=1.5, ",alpha,"=5")),expression(paste(m,"=2.0, ",alpha,"=50")),expression(paste(m,"=2.5, ",alpha,"=100")))) # Figure 6.4 # Page 7 # Left y = seq(-3.5,3.5,0.01) pdf = dnorm(y,0,1) plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxp=c(-3,3,6),cex.lab=1.25) abline(h=0) # Right u = seq(0,9,0.001) pdf = dchisq(u,1) plot(u,pdf,type="l",xlab="u",ylab="PDF",xlim=c(0,9),ylim=c(0,2),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) # Figure 6.5 # Page 8 x <- seq(0,2,0.001) y1 <- rep(2,length(x)) y2 <- rep(0,length(x)) plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(0,2),labels=c("0","")) axis(side=2,at=c(0,2),labels=c("0","")) # Figure 6.6 # Page 9 x <- seq(0,3,0.01) y1 <- 3-x y2 <- rep(0,length(x)) plot(x,y1,type="l",bty="n",xlim=c(0,4),ylim=c(0,4),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(0,3,5),labels=c("0","u","")) axis(side=2,at=c(0,3,5),labels=c("0","u","")) # Figure 6.7 # Page 12 # Left y = seq(0,1,0.001) pdf = dbeta(y,2,6) plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25) abline(h=0) # Right u = seq(0,1,0.001) pdf = dbeta(u,6,2) plot(u,pdf,type="l",xlab="u",ylab="PDF",cex.lab=1.25) abline(h=0) # Figure 6.8 # Page 13 # Left y = seq(-pi/2,pi/2,0.01) pdf = (1+y-y)/pi matplot(y,pdf,xlab="y",ylab="PDF",type="l",xaxt="n",yaxt="n",xlim=c(-2.5,2.5),ylim=c(0,0.75),cex.lab=1.25) axis(1, at = c(-pi/2, -pi/4, 0, pi/4, pi/2),labels = expression(-pi/2, -pi/4, 0, pi/4, pi/2)) axis(2, at = c(0, 1/pi),labels = expression(0, 1/pi)) lines(c(-2.5,-pi/2),c(0,0),lty=1) lines(c(pi/2,2.5),c(0,0),lty=1) # Right y = seq(-pi/2,pi/2,0.01) u = tan(y) matplot(y,u,xlab="y",ylab="h(y)",type="l",xaxt="n",yaxt="n",xlim=c(-2.5,2.5),ylim=c(-10,10),cex.lab=1.25) abline(v=0,lty=2) abline(v=-pi/2,lty=2) abline(v=pi/2,lty=2) axis(1, at = c(-pi/2, -pi/4, 0, pi/4, pi/2),labels = expression(-pi/2, -pi/4, 0, pi/4, pi/2)) axis(2, at = c(-10,-5,0,5,10),labels = c("-10","-5","0","5","10")) # Figure 6.9 # Page 14 u = seq(-8,8,0.001) pdf = (1/pi)*(1/(1+u^2)) plot(u,pdf,type="l",xlab="u",ylab="PDF",xaxp=c(-8,8,8),cex.lab=1.25) abline(h=0) # Figure 6.10 # Page 16 u = seq(0,30,0.01) plot(u,dlnorm(u,1,1),type="l",lty=1,xlab="u",ylab="PDF",cex.lab=1.25) lines(u,dlnorm(u,1.5,0.5),lty=4) lines(u,dlnorm(u,2.5,0.25),lty=8) abline(h=0) legend(15,0.20,lty=c(1,4,8),c(expression(paste(mu,"=1, ",sigma,"=1")), expression(paste(mu,"=1.5, ",sigma,"=0.5")),expression(paste(mu,"=2.5, ",sigma,"=0.25 ")))) # Figure 6.11 # Page 26 # Left x <- seq(0,2,0.001) y1 <- rep(2,length(x)) y2 <- rep(0,length(x)) plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(0,2),labels=c("0","")) axis(side=2,at=c(0,2),labels=c("0","")) # Right x <- seq(0,2,0.001) y1 <- rep(1/4,length(x)) y2 <- rep(0,length(x)) plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(u[1]),ylab=expression(u[2]),cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(0,2),labels=c("0","")) axis(side=2,at=c(0,1/4,2),labels=c("0","1","")) # Figure 6.12 # Page 29 u = seq(-6,6,0.001) pdf = 0.5*exp(-abs(u)) plot(u,pdf,type="l",xlab=expression(u[2]),ylab="PDF",ylim=c(0,0.5),cex.lab=1.25) abline(h=0) # Figure 6.13 # Page 30 # Left x <- seq(0,2,0.001) y1 <- rep(2,length(x)) y2 <- rep(0,length(x)) plot(x,y1,type="l",bty="n",xlim=c(0,2),ylim=c(0,2),xaxt='n',yaxt='n',xlab=expression(y[1]),ylab=expression(y[2]),cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(0,2),labels=c("0","")) axis(side=2,at=c(0,2),labels=c("0","")) # Right x <- seq(0,2,0.001) y1 <- x y2 <- -x plot(x,y1,type="l",bty="n",xlim=c(-2,2),ylim=c(-2,2),xaxt='n',yaxt='n',xlab=expression(u[1]),ylab=expression(u[2]),cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(-2,0,2),labels=c("","0","")) axis(side=2,at=c(-2,0,2),labels=c("","0","")) # Figure 6.14 # Page 38 y = seq(0,80,0.1) pdf = dexp(y,1/12) plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25) abline(h=0) abline(v=0,lty=2) # Figure 6.15 # Page 39 # Left y = seq(0,6,0.1) pdf = dexp(y,1/(12/14)) plot(y,pdf,type="l",xlab="y",ylab="PDF (minimum)",cex.lab=1.25) abline(h=0) abline(v=0,lty=2) # Right y = seq(0,120,0.1) pdf = (14/12)*exp(-y/12)*(1-exp(-y/12))^13 plot(y,pdf,type="l",xlab="y",ylab="PDF (maximum)",cex.lab=1.25) abline(h=0) # Figure 6.16 # Page 41 # Upper left y = seq(0,1,0.01) pdf = rep(1,length(y)) plot(y,pdf,type="l",xlab="y",ylab="PDF",yaxt='n',xlim=c(0,1),ylim=c(0,3/2),cex.lab=1.25) axis(side=2,at=c(0,1),labels=c("0","1")) abline(h=0) lines(c(1,1.5),c(0,0),lty=1) # Other figures n = 25 # number of observations k = 13 # order statistic position y = seq(0,1,0.001) plot(y,dbeta(y,k,n-k+1),type="l",lty=1,xlab="y",ylab="PDF (median)",ylim=c(0,max(dbeta(y,k,n-k+1))),cex.lab=1.25) abline(h=0) # abline(v=0,lty=2) # use if k=1 (min) # abline(v=1,lty=2) # use if k=n (max) # Figure 6.17 # Page 43 # Left x <- seq(0,1,0.001) y1 <- x y2 <- rep(1,length(x)) plot(x,y1,type="l",bty="n",xlim=c(0,1.1),ylim=c(0,1.1),xlab=expression(y[1]),ylab=expression(y[n]),xaxt='n',yaxt='n',cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(0,1,1.2),labels=c("0","1","")) axis(side=2,at=c(0,1,1.2),labels=c("0","1","")) # Right x <- seq(0,1,0.001) y1 <- x y2 <- rep(1,length(x)) plot(x,y1,type="l",bty="n",xlim=c(0,1.1),ylim=c(0,1.1),xlab=expression(u[1]),ylab=expression(u[2]),xaxt='n',yaxt='n',cex.lab=1.25) polygon(c(x,rev(x)),c(y2,rev(y1)),col="grey") axis(side=1,at=c(0,1,1.2),labels=c("0","1","")) axis(side=2,at=c(0,1,1.2),labels=c("0","1",""))