###############################################
## Author: Joshua M. Tebbs
## Date: 27 February 2019
## Update: 4 April 2019
## STAT 512 course notes: R Code
## Chapter 9
###############################################

# Figure 9.1
# Page 119
options("scipen"=0, "digits"=22)
accidents = c(rep(0,32),rep(1,26),rep(2,12),rep(3,7),rep(4,4),rep(5,2),rep(6,1))
sum(accidents)
prod(factorial(accidents))
theta = seq(0,3,0.01)
like = (theta^103)*exp(-84*theta)/prod(factorial(accidents))
plot(theta,like,type="l",lty=1,xlab=expression(theta),ylab="Likelihood function",yaxt="n",xaxp=c(0,3,3),cex.lab=1.25)
abline(h=0)


# Figure 9.2
# Page 121
bulbs = c(4.43,5.93,3.74,5.82,5.90,2.90,2.64,6.49,5.31,8.49,
	1.01,1.07,1.41,3.42,1.22,4.01,0.57,1.47,2.81,8.52,
	0.52,4.77,0.85,2.21,6.85,3.43,1.87,5.15,2.02,10.58)
sum(bulbs^2)
prod(bulbs)
theta = seq(0,50,0.01)
c = (2^30)*prod(bulbs)
like = c*(theta^(-30))*exp(-645.0/theta)
plot(theta,like,type="l",lty=1,xlab=expression(theta),ylab="Likelihood function",yaxt="n",xaxp=c(0,50,5),cex.lab=1.25)
abline(h=0)


# Figure 9.3
# Page 124
y = seq(2,13,0.1)
pdf = dexp(y-2,1/2)
plot(y,pdf,type="l",xlab="y",ylab="PDF",xlim=c(0,13),xaxt="n",yaxt="n",cex.lab=1.25)
abline(h=0)
abline(v=2,lty=2)
axis(side=1,at=c(0,2),labels=expression(0,theta))


# Figure 9.4
# Page 143
options("scipen"=0, "digits"=22)
accidents = c(rep(0,32),rep(1,26),rep(2,12),rep(3,7),rep(4,4),rep(5,2),rep(6,1))
sum(accidents)
theta = seq(0,3,0.01)
like = (theta^103)*exp(-84*theta)/prod(factorial(accidents))
plot(theta,like,type="l",lty=1,xlab=expression(theta),ylab="Likelihood function",yaxt="n",xaxp=c(0,3,3),cex.lab=1.25)
abline(h=0)
points(x=mean(accidents),y=0,pch=19,cex=1.25)


# Figure 9.5
# Page 145
bulbs = c(4.43,5.93,3.74,5.82,5.90,2.90,2.64,6.49,5.31,8.49,
	1.01,1.07,1.41,3.42,1.22,4.01,0.57,1.47,2.81,8.52,
	0.52,4.77,0.85,2.21,6.85,3.43,1.87,5.15,2.02,10.58)
sum(bulbs^2)
prod(bulbs)
theta = seq(0,50,0.01)
c = (2^30)*prod(bulbs)
like = c*(theta^(-30))*exp(-645.0/theta)
plot(theta,like,type="l",lty=1,xlab=expression(theta),ylab="Likelihood function",yaxt="n",xaxp=c(0,50,5),cex.lab=1.25)
abline(h=0)
points(x=sum(bulbs^2)/length(bulbs),y=0,pch=19,cex=1.25)


# Figure 9.6
# Page 146
theta = seq(1,3,0.01)
like = 1/theta^5
plot(theta,like,type="l",lty=1,xlab=expression(theta),ylab="Likelihood function",yaxt="n",xaxt="n",
xlim=c(0,3),ylim=c(0,1.1),cex.lab=1.25)
abline(h=0)
points(x=1,y=0,pch=19,cex=1.25)
axis(side=1,at=c(0,1),labels=expression(0,y[(n)]))


# Figure 9.7
# Page 150
max.rain = read.table("C:\\Users\\tebbs.MATHSTAT\\Documents\\texfiles\\Classes\\USC\\stat512\\s19\\data\\rainfall.txt",header=TRUE)
max.rain = max.rain[,1]
bins = seq(0,70,5)
hist(max.rain,prob=TRUE,breaks=bins,xlab="Maximum rainfall (inches)",ylab="Proportion",main="",col="lightblue",cex.lab=1.25)


# Example 9.20
# Rainfall data analysis
# Pages 150-152
n = 142 # sample size
m2 = (1/n)*sum(max.rain^2) # second sample moment; needed for MOM
# MOM estimates (see Example 9.16 notes)
alpha.mom = (mean(max.rain))^2/(m2-(mean(max.rain))^2)
beta.mom = (m2-(mean(max.rain))^2)/mean(max.rain)
# Sufficient statistics
t1 = sum(log(max.rain))
t2 = sum(max.rain)
# Negative loglikelihood function (to be minimised)
# x1 = alpha
# x2 = beta
loglike = function(x){
	x1<-x[1] 
	x2<-x[2]
	n*log(gamma(x1))+n*x1*log(x2)-t1*(x1-1)+t2/x2
	}
# Use "optim" function to maximise the loglikelihood function
mle = optim(par=c(alpha.mom,beta.mom),fn=loglike,method="Nelder-Mead")
c(alpha.mom,beta.mom) # MOM estimates
mle$par # MLEs


# Figure 9.8
# Page 152
# Left 
max.rain = read.table("C:\\Users\\tebbs.MATHSTAT\\Documents\\texfiles\\Classes\\USC\\stat512\\s19\\data\\rainfall.txt",header=TRUE)
max.rain = max.rain[,1]
bins = seq(0,70,5)
hist(max.rain,prob=TRUE,breaks=bins,xlab="Maximum rainfall (inches)",ylab="Proportion",main="",col="lightblue",cex.lab=1.25)
lines(sort(max.rain),dgamma(sort(max.rain),mle$par[1],1/mle$par[2]),col="red",lwd=2)
# Right 
plot(qgamma(ppoints(max.rain),mle$par[1],1/mle$par[2]),sort(max.rain),pch=16,xlab="Gamma percentiles",ylab="Observed data",col="black")
abline(a=0,b=1,lty=2)