### MME fitted gamma

# X_1, X_2, ...,X_n iid Gamma(shape=,rate=lambda), find the MMEs.
# p.50 f(x|alpha,lambda)= lambda^alpha/gamma(alpha)) x^(alpha-1) e^-(lambda*x)

# parameters

n <- 100

alpha <- .5; a.min <- 0.1; a.max <- 2;
lambda <- 2; l.max <- 10


# generate random data

x <- rgamma(n, shape = alpha, rate = lambda)

x.orig <- x	# save a copy of the original data

# relative frequency histogram of the simulated data

hist(x,probability=T)

x.max <- max(hist(x)$breaks)
y.max <- max(hist(x)$density)

# since X is continuous, really should use a density estimator

plot(density(x, n=50, window="g"),type="l",xlab="x",ylab="density")

# MMEs p.249

sigma2.hat <- (1/n)*sum( (x - mean(x))^2 )

alpha.hat.mme <- mean(x)^2/sigma2.hat
lambda.hat.mme <- mean(x)/sigma2.hat

# plot the fitted model to the simulated data

xx <- seq(0,x.max,.1)

yy <- dgamma(xx, shape = alpha.hat.mme, rate = lambda.hat.mme)

plot(xx,yy, main="MME fitted gamma",type="l")

# make plot of histogram with the fitted model

hist(x,probability=T,xlim=c(0,x.max),ylim=c(0,y.max),main="MME fitted gamma",
  xlab="x",ylab="density")
par(new=T)
plot(xx,yy, type="l",xlim=c(0,x.max),ylim=c(0,y.max),xlab="", ylab="")

# make a plot of the density estimator with the fitted model.

plot(density(x, n=50, window="g"),type="l",xlim=c(0,x.max),ylim=c(0,y.max),
   main="MME fitted gamma",xlab="x",ylab="density")
par(new=T)
plot(xx,yy, type="l",xlim=c(0,x.max),ylim=c(0,y.max),
   xlab="x",ylab="density")

# parametric bootstrap the MME

B <- 1000

alpha.hat.mme.star <- numeric(B)
lambda.hat.mme.star <- numeric(B)

for(i in 1:B){
	x <- rgamma(n,shape=alpha.hat.mme,rate=lambda.hat.mme)
	sigma2.hat.star <- (1/n)*sum( (x - mean(x))^2 )
	alpha.hat.mme.star[i] <- mean(x)^2/sigma2.hat.star
	lambda.hat.mme.star[i] <- mean(x)/sigma2.hat.star
}

plot(density(alpha.hat.mme.star, n=50, window="g"),type="l",xlab="x",ylab="density")

plot(density(lambda.hat.mme.star, n=50, window="g"),type="l",xlab="x",ylab="density")

alpha.mme.boot.m <- mean(alpha.hat.mme.star)
alpha.mme.boot.s <- sd(alpha.hat.mme.star)
alpha.mme.boot.ci <- quantile(alpha.hat.mme.star,c(0.025,0.975))

lambda.mme.boot.m <- mean(lambda.hat.mme.star)
lambda.mme.boot.s <- sd(lambda.hat.mme.star)
lambda.mme.boot.ci <- quantile(lambda.hat.mme.star,c(0.025,0.975))

#############################################################################

### MLE fitted gamma

# Note:  There is a function is R called the digamma that will calculate
# the derivative of the log of the gamma function.  

# digamma(x) = psi(x) = d/dx {ln Gamma(x)} = Gamma'(x) / Gamma(x)

alpha.d <- function(a){
	 n*log(a) - n*log(mean(x)) + sum(log(x)) - n*digamma(a)
}

# find the zero of the derivative

# direct method, plot the function and find where it crosses zero
# try different min and max values to locate the zero

a.min 
a.max 

a <- seq(a.min,a.max,0.01)

pp <- n*log(a) - n*log(mean(x)) + sum(log(x)) - n*digamma(a)

A <- matrix(c(a,pp),ncol=2)
plot(A,type="l",xlab="alpha",ylab="alpha.d")

X <- matrix(NA,nrow=length(a),ncol=2)
X[,1] <- a
X[,2] <- 0

lines(X)

# Find the value of a where the derivative function is 
# just above zero

A1 <- matrix(A[A[,2]>0],ncol=2)
alpha.mle1 <- A1[length(A1[,1]),1]

# Find the value of a where the derivative function is
# just below zero

A2 <- matrix(A[A[,2]<0],ncol=2)
alpha.mle2 <- A2[1,1]

# Take the mle to be the interpolated value

alpha.mle <- (alpha.mle1 + alpha.mle2)/2
alpha.mle

lambda.mle <- alpha.mle/mean(x)
lambda.mle

# using the uniroot function in R to find the zero of the derivative

alpha.hat.mle <- uniroot(alpha.d, c(a.min,a.max))

alpha.hat.mle <- alpha.hat.mle$root
alpha.hat.mle

lambda.hat.mle <- alpha.hat.mle /mean(x)
lambda.hat.mle

# plot the fitted model to the simulated data

xx <- seq(0,x.max,.1)

yy <- dgamma(xx, shape = alpha.hat.mle, rate = lambda.hat.mle)

plot(xx,yy, main="MLE fitted gamma",type="l")

# make plot of histogram with the fitted model

hist(x,probability=T,xlim=c(0,x.max),ylim=c(0,y.max),main="MLE fitted gamma",
  xlab="x",ylab="density")
par(new=T)
plot(xx,yy, type="l",xlim=c(0,x.max),ylim=c(0,y.max),xlab="", ylab="")

# make a plot of the density estimator with the fitted model.

plot(density(x, n=50, window="g"),type="l",xlim=c(0,x.max),ylim=c(0,y.max),
   main="MLE fitted gamma",xlab="x",ylab="density")
par(new=T)
plot(xx,yy, type="l",xlim=c(0,x.max),ylim=c(0,y.max),
   xlab="x",ylab="density")

# parametric bootstrap the MLE

B <- 1000

alpha.hat.mle.star <- numeric(B)
lambda.hat.mle.star <- numeric(B)

a.min <- 0.1
a.max <- 2

for(i in 1:B){
        x <- rgamma(n,shape=alpha.hat.mle,rate=lambda.hat.mle)
        sigma2.hat.star <- (1/n)*sum( (x - mean(x))^2 )
        alpha.hat.mle.star[i] <- uniroot(alpha.d, c(a.min,a.max))$root
        lambda.hat.mle.star[i] <- mean(x)/sigma2.hat.star
}

plot(density(alpha.hat.mle.star, n=50, window="g"),type="l",xlab="x",ylab="density")

plot(density(lambda.hat.mle.star, n=50, window="g"),type="l",xlab="x",ylab="density")

alpha.mle.boot.m <- mean(alpha.hat.mle.star)
alpha.mle.boot.s <- sd(alpha.hat.mle.star)
alpha.mle.boot.ci <- quantile(alpha.hat.mle.star,c(0.025,0.975))

lambda.mle.boot.m <- mean(lambda.hat.mle.star)
lambda.mle.boot.s <- sd(lambda.hat.mle.star)
lambda.mle.boot.ci <- quantile(lambda.hat.mle.star,c(0.025,0.975))

#############################################################################
# make and print tables for the output

MME <- c(alpha.hat.mme,alpha.mme.boot.s,lambda.hat.mme,lambda.mme.boot.s)
MLE <- c(alpha.hat.mle,alpha.mle.boot.s,lambda.hat.mle,lambda.mle.boot.s)

rlabel <- c("alpha.hat","alpha.hat.boot.se","lambda.hat","lambda.hat.boot.se")
clabel <- c("mme", "mle")

dm <- list(rlabel, clabel)

M1 <- matrix(c(MME,MLE),ncol=2,dimnames=dm)


MME <- c(alpha.mme.boot.ci,lambda.mme.boot.ci)
MLE <- c(alpha.mle.boot.ci,lambda.mle.boot.ci)

rlabel <- c("alpha.hat","lambda.hat")
clabel <- c("mme 0.025","mme 0.975","mle 0.025","mle 0.975")

dm <- list(rlabel, clabel)

M2 <- matrix(c(MME,MLE),ncol=4,dimnames=dm)

MME <- c(diff(alpha.mme.boot.ci),diff(lambda.mme.boot.ci))
MLE <- c(diff(alpha.mle.boot.ci),diff(lambda.mle.boot.ci))

rlabel <- c("alpha.hat","lambda.hat")
clabel <- c("mme ci length","mle ci length")

dm <- list(rlabel, clabel)

M3 <- matrix(c(MME,MLE),ncol=2,dimnames=dm)

M1
M2
M3

#######################################################################################

# the 3d likelihood for the gamma

par(mfrow=c(2,2))

alp <- seq(0.1,a.max,0.05)
lam <- seq(0.1,l.max,0.05)

log.lik.gam <- function(alp,lam){
	n*alp*log(lam) + (alp - 1)*sum(log(x.orig)) - lam*sum(x.orig) - n*log(gamma(alp))
}


zz <- outer(alp,lam,log.lik.gam)

persp(alp, lam, zz,theta = 30, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "alpha", ylab = "lambda", zlab = "likelihood")


zz <- log.lik.gam(alpha.hat.mle,lam)
plot(lam,zz,type="l",main="for alpha.mle",xlab="lambda",ylab="lik")

zz <- log.lik.gam(alp,lambda.hat.mle)
plot(alp,zz,type="l",main="for lambda.mle",xlab="alpha",ylab="lik")