#Find gamma MLEs options(width = 60) ################################################################################ # Simulate data set.seed(4869) n <- 30 alpha <- 7/2 beta <- 2 #Chisq(7) y <- rgamma(n = n, shape = alpha, scale = beta) round(y,2) win.graph(width = 12, pointsize = 16) par(mfrow = c(1,3)) #PDF hist(x = y, xlab = "x", freq = FALSE, ylim = c(0, 0.13), main = "Histogram", xlim = c(0,22)) curve(expr = dgamma(x = x, shape = alpha, scale = beta), col = "red", add = TRUE) #CDFs plot.ecdf(x = y, verticals = TRUE, do.p = FALSE, main = "CDFs", lwd = 2, col = "black", xlab = "x", ylab = "CDF") curve(expr = pgamma(q = x, shape = alpha, scale = beta), col = "red", add = TRUE, lwd = 2, lty = "dotted") # QQ-Plot gamma.quant <- qgamma(p = seq(from = 1/(n+1), to = 1-1/(n+1), by = 1/(n+1)), shape = alpha, scale = beta) plot(x = sort(y), y = gamma.quant, main = "QQ-plot", xlab = "x", ylab = "Gamma quantiles") abline(a = 0, b = 1) ################################################################################ # Maximum likelihood estimation logL1 <- function(theta, y) { n <- length(y) alpha <- theta[1] beta <- theta[2] sum(dgamma(x = y, shape = alpha, scale = beta, log = TRUE)) } alpha.mom <- mean(y)^2/((n-1)*var(y)/n) beta.mom <- ((n-1)*var(y)/n)/mean(y) data.frame(alpha.mom, beta.mom) gam.mle1 <- optim(par = c(alpha.mom, beta.mom), fn = logL1, y = y, method = "BFGS", control = list(fnscale = -1), hessian = TRUE) gam.mle1 #Estimated covariance matrix cov.mat <- -solve(gam.mle1$hessian) cov.mat data.frame(name = c("alpha", "beta"), estimate = gam.mle1$par, SE = sqrt(diag(cov.mat))) ################################################################################ # Contour plot #Set up sequence of parameter values# alpha<-seq(0.25,10,0.1) beta<-seq(0.25,10,0.1) pars<-as.matrix(expand.grid(alpha, beta)) #Evaluate logL values# loglik<-matrix(data = NA, nrow = nrow(pars), ncol = 1) for(i in 1:nrow(pars)){ theta<-pars[i,] loglik[i,1]<-c(logL1(theta = theta, y = y)) } #Put log(L) into matrix - need to be careful about ordering loglik.mat<-matrix(data = loglik, nrow = length(alpha), ncol = length(beta), byrow = FALSE) head(loglik) loglik.mat[1:5, 1:5] logL1(theta = c(0.25, 0.25), y = y) logL1(theta = c(0.25, 0.35), y = y) logL1(theta = c(0.35, 0.25), y = y) win.graph() par(pty = "s", mfrow = c(1,1)) contour(x = alpha, y = beta, z = loglik.mat, levels = c(-84, -85, -86, -88, -90, -100, -200, -300, -400), xlab = expression(alpha), ylab = expression(beta), main = "Contour plot of log(L)") abline(h = gam.mle1$par[2], col = "blue") abline(v = gam.mle1$par[1], col = "blue") max(loglik) #From help for contour: Notice that contour interprets the z matrix as a table of f(x[i], y[j]) values, # so that the x axis corresponds to row number and the y axis to column number, with column 1 at the # bottom, i.e. a 90 degree counter-clockwise rotation of the conventional textual layout. library(package = rgl) #Package that does 3D interactive plots open3d() #Open plot window #3D plot with gridlines persp3d(x = alpha, y = beta, z = loglik.mat, xlab = "alpha", ylab = "beta", zlab = "log(L)", ticktype = "detailed", col = "red") spheres3d(x = gam.mle1$par[1], y = gam.mle1$par[2], z = max(loglik), radius = 10) lines3d(x = c(gam.mle1$par[1], gam.mle1$par[1]), y = c(gam.mle1$par[2], gam.mle1$par[2]), z = c(min(loglik), max(loglik))) grid3d(c("x", "y+", "z")) #Plot with different scale alpha<-seq(1,3.5,0.1) beta<-seq(1,3.5,0.1) pars<-as.matrix(expand.grid(alpha, beta)) loglik<-matrix(data = NA, nrow = nrow(pars), ncol = 1) for(i in 1:nrow(pars)){ theta<-pars[i,] loglik[i,1]<-c(logL1(theta = theta, y = y)) } loglik.mat<-matrix(data = loglik, nrow = length(alpha), ncol = length(beta), byrow = FALSE) persp3d(x = alpha, y = beta, z = loglik.mat, xlab = "alpha", ylab = "beta", zlab = "log(L)", ticktype = "detailed", col = "red") spheres3d(x = gam.mle1$par[1], y = gam.mle1$par[2], z = max(loglik), radius = 5) lines3d(x = c(gam.mle1$par[1], gam.mle1$par[1]), y = c(gam.mle1$par[2], gam.mle1$par[2]), z = c(min(loglik), max(loglik))) grid3d(c("x", "y+", "z")) ################################################################################ # Maximum likelihood estimation using optim() with L-BFGS-B gam.mle1.constr <- optim(par = c(alpha.mom, beta.mom), fn = logL1, y = y, method = "L-BFGS-B", control = list(fnscale = -1), hessian = TRUE, lower = c(0, 0), upper = c(Inf, Inf)) gam.mle1.constr ################################################################################ # Maximum likelihood estimation using log(par) logL2 <- function(theta, y) { n <- length(y) alpha <- exp(theta[1]) beta <- exp(theta[2]) sum(dgamma(x = y, shape = alpha, scale = beta, log = TRUE)) } gam.mle2 <- optim(par = log(c(alpha.mom, beta.mom)), fn = logL2, y = y, method = "BFGS", control = list(fnscale = -1), hessian = TRUE) gam.mle2 exp(gam.mle2$par) #Estimated covariance matrix cov.mat.log <- -solve(gam.mle2$hessian) cov.mat.log gdot <- matrix(data = c(exp(gam.mle2$par[1]), 0, 0, exp(gam.mle2$par[2])), nrow = 2, ncol = 2) cov.mat2 <- gdot %*% cov.mat.log %*% gdot cov.mat2 data.frame(name = c("alpha", "beta"), estimate = exp(gam.mle2$par), SE = sqrt(diag(cov.mat2))) ################################################################################ # Maximum likelihood estimation using poor coding logL3 <- function(theta, y) { n <- length(y) alpha <- theta[1] beta <- theta[2] -n*log(gamma(alpha)) - n*alpha*log(beta) - sum(y)/beta + (alpha-1)*sum(log(y)) } #Notice tried different starting values gam.mle3 <- optim(par = c(1, 1), fn = logL3, y = y, method = "BFGS", control = list(fnscale = -1), hessian = TRUE) #gam.mle3 #Estimated covariance matrix cov.mat3 <- -solve(gam.mle3$hessian) cov.mat3 data.frame(name = c("alpha", "beta"), estimate = gam.mle3$par, SE = sqrt(diag(cov.mat3))) #Problems do not occur with logL2() gam.mle4 <- optim(par = log(c(1, 1)), fn = logL2, y = y, method = "BFGS", control = list(fnscale = -1), hessian = TRUE) cov.mat.log4 <- -solve(gam.mle4$hessian) gdot4 <- matrix(data = c(exp(gam.mle4$par[1]), 0, 0, exp(gam.mle4$par[2])), nrow = 2, ncol = 2) cov.mat4 <- gdot4 %*% cov.mat.log4 %*% gdot4 cov.mat4 data.frame(name = c("alpha", "beta"), estimate = exp(gam.mle4$par), SE = sqrt(diag(cov.mat4))) #