############################################### ## 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)