############################################### ## Author: Joshua M. Tebbs ## Date: 21 August 2018 ## Update: 10 November 2018 ## STAT 511 course notes: R Code ## Chapter 4 ############################################### # Figure 4.1 # Page 74 # Left y = c(0,1,2,3) prob = dbinom(y,3,0.4) plot(y,prob,type="h",xlab="y",ylab="PMF",ylim=c(0,max(prob)),xaxt='n',cex.lab=1.25) abline(h=0) axis(side=1,at=y,labels=c("0","1","2","3")) # Right y = c(0,1,2,3) cdf = c(0,cumsum(prob)) cdf.plot = stepfun(y,cdf,f=0) plot.stepfun(cdf.plot,xlab="y",ylab="CDF",verticals=FALSE,do.points=TRUE,main="",pch=16,cex.lab=1.25) # Figure 4.2 # Page 77 y = seq(0,4,0.01) cdf = 1-exp(-y^2) plot(y,cdf,type="l",xlab="y",ylab="CDF",xlim=c(-1,4),ylim=c(0,1),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) # Figure 4.3 # Page 78 bmi = read.table("C:\\Users\\Tebbs\\Documents\\texfiles\\Classes\\USC\\stat511\\f18\\data\\bmi.txt",header=TRUE) bmi = bmi[,1] bins = seq(10,50,2) hist(bmi,prob=TRUE,breaks=bins,xlab="Body mass index (BMI)",ylab="Proportion",main="",col="lightblue",cex.lab=1.25) lines(density(bmi),col="red",lwd=2) # Figure 4.4 # Page 80 # Left y = seq(0,2,0.01) pdf = (3/8)*(4*y-2*y^2) plot(y,pdf,type="l",xlab="y",ylab="PDF",xlim=c(-1,3),cex.lab=1.25) abline(h=0) # Right y = seq(0,2,0.01) cdf = c((3/8)*(2*y^2 - (2/3)*y^3)) plot(y,cdf,type="l",xlab="y",ylab="CDF",xlim=c(-1,3),ylim=c(0,1),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) lines(c(2,3),c(1,1),lty=1) # Figure 4.5 # Page 81 # Left y = seq(0,2,0.01) pdf = (3/8)*(4*y-2*y^2) plot(y,pdf,type="l",xlab="y",ylab="PDF",xlim=c(-1,3),cex.lab=1.25) abline(h=0) x = seq(0,0.5,0.001) y = (3/8)*(4*x-2*x^2) polygon(c(0,x,0.5),c(0,y,0),col="lightblue") points(x=0.5,y=0,pch=19,cex=1.5) # Right y = seq(0,2,0.01) cdf = c((3/8)*(2*y^2 - (2/3)*y^3)) plot(y,cdf,type="l",xlab="y",ylab="CDF",xlim=c(-1,3),ylim=c(0,1),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) lines(c(2,3),c(1,1),lty=1) points(x=0.5,y=(3/8)*(2*(0.5)^2-(2/3)*(0.5)^3),pch=19,cex=1.5) # Figure 4.6 # Page 82 # Left y = seq(-7,7,0.01) pdf = exp(-y)/(1+exp(-y))^2 plot(y,pdf,type="l",xlab="y",ylab="PDF",xlim=c(-7,7),cex.lab=1.25) abline(h=0) x = seq(-2,2,0.001) y = exp(-x)/(1+exp(-x))^2 polygon(c(-2,x,2),c(0,y,0),col="lightblue") points(x=-2,y=0,pch=19,cex=1.5) points(x=2,y=0,pch=19,cex=1.5) # Right y = seq(-7,7,0.01) cdf = 1/(1+exp(-y)) plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25) abline(h=0) points(x=-2,y=1/(1+exp(2)),pch=19,cex=1.5) points(x=2,y=1/(1+exp(-2)),pch=19,cex=1.5) # Figure 4.7 # Page 84 # Left y = seq(0,35,0.01) pdf = 3000/(10+y)^4 plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25) abline(h=0) abline(v=0,lty=2) x = seq(0,10/((0.5)^(1/3))-10,0.001) y = 3000/(10+x)^4 polygon(c(0,x,10/((0.5)^(1/3))-10),c(0,y,0),col="lightblue") points(x=10/((0.5)^(1/3))-10,y=0,pch=19,cex=1.5) # Right y = seq(0,35,0.01) cdf = 1-(10/(10+y))^3 plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) abline(h=0.5,lty=2) points(x=10/((0.5)^(1/3))-10,y=0,pch=19,cex=1.5) # Figure 4.8 # Page 86 y = seq(0,10,0.01) pdf = y*exp(-y) plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25) abline(h=0) points(x=2,y=0,pch=19,cex=1.5) # Figure 4.9 # Page 89 y = seq(0,6,0.01) plot(y,dexp(y,1/0.5),type="l",lty=4,xlab="y",ylab="PDF",cex.lab=1.25) lines(y,dexp(y,1),lty=1) lines(y,dexp(y,1/2),lty=8) abline(h=0) abline(v=0,lty=2) legend(3,1,lty = c(4,1,8),c(expression(paste(beta, "=0.5")),expression(paste(beta, "=1")),expression(paste(beta, "=2")))) # Figure 4.10 # Page 91 # Left y = seq(0,2*pi,0.01) pdf = (1+y-y)/(2*pi) matplot(y,pdf,xlab="y",ylab="PDF",type="l",xaxt="n",yaxt="n",xlim=c(0,7.5),ylim=c(0,0.5),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) axis(1, at = c(0, pi/2, pi, 3*pi/2, 2*pi),labels = expression(0, pi/2, pi, 3*pi/2, 2*pi)) axis(2, at = c(0, 1/(2*pi)),labels = expression(0, 1/2*pi)) # Right cdf = y/(2*pi) matplot(y,cdf,xlab="y",ylab="CDF",type="l",xaxt="n",xlim=c(0,7.5),ylim=c(0,1),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) lines(c(2*pi,7.5),c(1,1),lty=1) axis(1, at = c(0, pi/2, pi, 3*pi/2, 2*pi),labels = expression(0, pi/2, pi, 3*pi/2, 2*pi)) # Figure 4.11 # Page 92 y = seq(-10,10,0.01) plot(y,dnorm(y,0,1),type="l",lty=1,xlab="y",ylab="PDF",cex.lab=1.25) lines(y,dnorm(y,-2,2),lty=4) lines(y,dnorm(y,1,3),lty=8) abline(h=0) legend(3.5,0.30,lty=c(1,4,8),c(expression(paste(mu,"=0, ",sigma,"=1")),expression(paste(mu,"=-2, ",sigma,"=2 ")),expression(paste(mu,"=1, ",sigma,"=3")))) # Figure 4.12 # Page 96 # Left y = seq(70,180,0.1) pdf = dnorm(y,125,15) plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxp=c(80,170,6),cex.lab=1.25) abline(h=0) x = seq(70,90,0.1) y = dnorm(x,125,15) polygon(c(70,x,90),c(0,y,0),col="lightblue") points(x=90,y=0,pch=19,cex=1.5) # Right y = seq(70,180,0.1) pdf = dnorm(y,125,15) plot(y,pdf,type="l",xlab="y",ylab="PDF",xaxp=c(80,170,6),cex.lab=1.25) abline(h=0) points(x=137.6243,y=0,pch=19,cex=1.5) # Figure 4.13 # Page 98 # Left y = seq(0,1200,0.1) pdf = dexp(y,1/190) plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25) abline(h=0) abline(v=0,lty=2) x = seq(100,1200,0.1) y = dexp(x,1/190) polygon(c(100,x,1200),c(0,y,0),col="lightblue") points(x=100,y=0,pch=19,cex=1.5) # Right y = seq(0,1200,0.1) cdf = pexp(y,1/190) plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25) abline(h=0) abline(v=0,lty=2) abline(h=0.9,lty=2) points(x=437.5,y=0,pch=19,cex=1.5) # Figure 4.14 # Page 101 y = seq(0,25,0.01) plot(y,dgamma(y,1.5,1/1.6),type="l",lty=1,xlab="y",ylab="PDF",cex.lab=1.25) lines(y,dgamma(y,2,1/2),lty=4) lines(y,dgamma(y,2.5,1/3),lty=8) abline(h=0) legend(10,0.20,lty=c(1,4,8),c(expression(paste(alpha,"=1.5, ",beta,"=1.6")),expression(paste(alpha,"=2.0, ",beta,"=2.0")),expression(paste(alpha,"=2.5, ",beta,"=3.0")))) # Figure 4.15 # Page 103 y = seq(0,1200,0.1) pdf = dgamma(y,2.7,1/100) plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25) abline(h=0) x = seq(365,1200,0.1) y = dgamma(x,2.7,1/100) polygon(c(365,x,1200),c(0,y,0),col="lightblue") points(x=365,y=0,pch=19,cex=1.5) # Figure 4.16 # Page 106 y = seq(0,1,0.01) plot(y,dbeta(y,3,3),type="l",lty=1,xlab="y",ylab="PDF",ylim=c(0,4),cex.lab=1.25) lines(y,dbeta(y,2,6.5),lty=4) lines(y,dbeta(y,4,1),lty=8) abline(h=0) abline(v=1,lty=2) legend(0.4,3.5,lty=c(1,4,8),c(expression(paste(alpha,"=3, ",beta,"=3")),expression(paste(alpha,"=2, ",beta,"=6.5")),expression(paste(alpha,"=4, ",beta,"=1")))) # Figure 4.17 # Page 108 y = seq(0,1,0.001) pdf = dbeta(y,1,5) plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25) abline(h=0) abline(v=0,lty=2) x = seq(qbeta(0.99,1,5),1,0.001) y = dbeta(x,1,5) polygon(c(qbeta(0.99,1,5),x,1),c(0,y,0),col="lightblue") text(0.7,0.3,0.01,cex = 1)