# PMD code for HW 1 # Problem 1: source('lnlmodels.r') hareM0.data <- c(6,66,14+26+18+24+21+30) hareM0.mle <- optim(c(70,0.9), lnlM0, data=hareM0.data, method='BFGS', control=list(fnscale=-1), hessian=T) hareMt.data <- c(6,66,14,26,18,24,21,30) hareMt.mle <- optim(c(79,rep(0.5,6)), lnlMt, data=hareMt.data, method='BFGS', control=list(fnscale=-1), hessian=T) hareMb.data <- c(6,66, 0+14+38+47+56+62, 0+2+9+15+15+26) hareMb.mle <- optim(c(75,0.5,0.3), lnlMb, data=hareMb.data, method='BFGS', control=list(fnscale=-1), hessian=T) hare.lnl <- sapply(list(hareM0.mle, hareMt.mle, hareMb.mle), function(x){x$value}) # extract the lnL value from each output list hare.AIC <- -2*hare.lnl + 2*c(2,7,3) Nhat <- 66:90 # potential range of ests of N hare.prof <- rep(NA, length(Nhat)) for (i in 1:length(Nhat)) { hare.prof[i] <- lnlMtp(Nhat[i], hareMt.data) } par(mar=c(5,4,0,0)+0.2) # shrink the upper and right margins plot(Nhat, hare.prof) Nhat[which.max(hare.prof)] # check that get same value for the mle from the profile function max(hare.prof) max(hare.prof)-qchisq(0.95,1)/2 abline(h=max(hare.prof)-qchisq(0.95,1)/2 ) which(hare.prof < (max(hare.prof)-qchisq(0.95,1)/2) ) # identify observations outside the ci, # want largest at the lower end and smallest at the upper end Nhat[c(3,19)] hare.est <- sapply(list(hareM0.mle, hareMt.mle, hareMb.mle), function(x){x$par[1]}) # extract population size from each output list # it is the first parameter dAIC <- hare.AIC - min(hare.AIC) wi <- exp(-dAIC/2) wi <- wi/sum(wi) hare.MAest <- sum(wi*hare.est) hare.MAest hare.vN <- sapply(list(hareM0.mle, hareMt.mle, hareMb.mle), function(x){ temp <- solve(-x\$hessian) temp[1,1] }) sum(wi*sqrt(hare.vN + (hare.est-hare.MAest)^2)) # problem 2: simulation code simN <- function(nsim, Nw, n) { # simulate nsim values of n_w and return vector of \hat{N}_{black} # using 'less biased' estimator for N_{black} phat <- Nw/(Nw+141) nw <- rbinom(10000, n, phat) (n+1)*(Nw+1)/(nw+1) - Nw - 1 } simN2 <- function(nsim, Nw, n) { # simulate nsim values of n_w and return vector of \hat{N}_{black} # using 'less biased' estimator for N_{black} phat <- Nw/(Nw+141) nw <- rbinom(10000, n, phat) nb <- n - nw (nb+1)*(Nw+1)/(nw+1) - 1 } sA <- simN(10000, 10, 20) sB <- simN(10000, 25, 20) sC <- simN(10000, 20, 10) # sA <- simN2(10000, 10, 20) # sB <- simN2(10000, 25, 20) # sC <- simN2(10000, 20, 10) # Note: variance of values from (2) is not exactly the same as # that from (1) because of the +1 terms. cat("Design A", mean(sA)-141, var(sA), sd(sA), '\n') cat("Design B", mean(sB)-141, var(sB), sd(sB), '\n') cat("Design C", mean(sC)-141, var(sC), sd(sC), '\n')