###############################################
## Author: Joshua M. Tebbs
## Date: 11 December 2011
## Update: 29 December 2017
## STAT 509 course notes: R Code
## Chapter 4
###############################################

# Example 4.1
# Page 40-41
# Continuous PDF and CDF
# Plot PDF
y = seq(0,1,0.01)
pdf = 3*y^2
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
abline(v=1,lty=2)
# Plot CDF
y0 = seq(-0.25,1.25,0.01)
y = seq(0,1,0.01)
cdf = c(rep(0,25),y^3,rep(1,25))
plot(y0,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)


# Example 4.2
# Page 43-44
# Continuous PDF and CDF
# Plot PDF
y = seq(12.5,13,0.01)
pdf = 20*exp(-20*(y-12.5))
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
abline(v=12.5,lty=2)
# Plot CDF
y = seq(12.5,13,0.01)
cdf = 1-exp(-20*(y-12.5))
plot(y,cdf,type="l",xlab="y",ylab="CDF",xlim=c(12.4,13),ylim=c(0,1),cex.lab=1.25)
abline(h=0)

# Calculate E(Y)
integrand.1 <- function(y){y*20*exp(-20*(y-12.5))}
integrate(integrand.1,lower=12.5,upper=Inf)
# Calculate var(Y)
integrand.2 <- function(y){(y-12.55)^2*20*exp(-20*(y-12.5))}
integrate(integrand.2,lower=12.5,upper=Inf)


# Figure 4.3
# Page 45
# Exponential PDFs
y = seq(0,8,0.01)
# Plot exponential pdf with lambda = 2
plot(y,dexp(y,2),type="l",lty=1,xlab="y",ylab="PDF")
# Add other pdfs
lines(y,dexp(y,1),lty=4)
lines(y,dexp(y,1/2),lty=8)
abline(h=0)
abline(v=0,lty=2)
# Add legend
legend(3,1,lty = c(1,4,8), 
	c(
	expression(paste(lambda, "=2")),
	expression(paste(lambda, "=1")),
	expression(paste(lambda, "=1/2"))
	))


# Example 4.3
# Page 46-47
# Exponential PDF and CDF
# Plot PDF
y = seq(0,15,0.01)
pdf = dexp(y,0.4)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
abline(v=0,lty=2)
# Plot CDF
cdf = pexp(y,0.4)
plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)


# Figure 4.5
# Page 50
# Gamma PDFs
y = seq(0,25,0.01)
# Plot gamma pdf with alpha = 1.5 and lambda = 0.6
plot(y,dgamma(y,1.5,0.6),type="l",lty=1,xlab="y",ylab="PDF")
# Add other pdfs
lines(y,dgamma(y,2,0.5),lty=4)
lines(y,dgamma(y,2.5,0.4),lty=8)
abline(h=0)
# Add legend
legend(10,0.20,lty = c(1,4,8), 
	c(
	expression(paste(alpha, "=1.5, ", lambda, "=0.6")),
	expression(paste(alpha, "=2.0, ", lambda, "=0.5")),
	expression(paste(alpha, "=2.5, ", lambda, "=0.4"))
	))


# Example 4.5
# Page 50-51
# Gamma PDF and CDF
# Plot PDF
y = seq(0,80,0.01)
pdf = dgamma(y,4,1/6)
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
# Plot CDF
cdf = pgamma(y,4,1/6)
plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)


# Figure 4.7
# Page 53
# Normal PDFs
y = seq(-10,10,0.01)
# Plot normal pdf with mu = 0 and sigma = 1 (standard normal pdf)
plot(y,dnorm(y,0,1),type="l",lty=1,xlab="y",ylab="PDF")
# Add other pdfs
lines(y,dnorm(y,-2,2),lty=4)
lines(y,dnorm(y,1,3),lty=8)
abline(h=0)
# Add legend
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"))
	))


# Example 4.6
# Page 53-54
# Normal PDF and CDF
# Plot PDF
y = seq(0,3,0.01)
pdf = dnorm(y,1.5,sqrt(0.16))
plot(y,pdf,type="l",xlab="y",ylab="PDF",cex.lab=1.25)
abline(h=0)
# Plot CDF
cdf = pnorm(y,1.5,sqrt(0.16))
plot(y,cdf,type="l",xlab="y",ylab="CDF",ylim=c(0,1),cex.lab=1.25)
abline(h=0)