# The purpose of this example is to make using nlm a bit easier.
# (1) We do not have to provide analytic derivatives.  It will
#     approximate them by finite differences if we don't.  This
#     won't work quite as well, but it will (usually) work.
# (2) We do not have to provide the data as an argument to mlogl.
#     We can just look them up in the R global environment.  Many
#     programmers scorn global variables and this should not be
#     done if you are writing code to be used by others.  But it
#     is OK for one-off projects.

# these are the data that were simulated in likelihood.R
nA <- 71
nB <- 23
nO <- 100
nAB <- 6

# this is the minimal mlogl that works
# Welllllll.  Taking out the stopifnot's would also work, but they
# are there for your protection.  If you leave them out you will
# probably use the function wrong an get errors that these checks
# would catch

mlogl <- function(theta) {
    stopifnot(is.numeric(theta))
    stopifnot(is.finite(theta))
    stopifnot(length(theta) == 2)
    alpha <- theta[1]
    beta <- theta[2]
    gamma <- 1 - alpha - beta
    if (alpha <= 0 | beta <= 0 | gamma <= 0) return(Inf)
    return(- (nA * log(alpha^2 + 2 * alpha * gamma) +
        nB * log(beta^2 + 2 * beta * gamma) +
        nO * log(gamma^2) +
        nAB * log(2 * alpha * beta)))

}

# everything from here on is just the same as in likelihood.R

nout <- nlm(mlogl, c(1/3, 1/3), hessian = TRUE)
nout <- nlm(mlogl, nout$estimate, hessian = TRUE)
nout$estimate

# was (0.21597208 0.07540475) in likelihood.Rout

alpha.hat <- nout$estimate[1]
beta.hat <- nout$estimate[2]
gamma.hat <- 1 - alpha.hat - beta.hat

# Wald type standard errors from inverse observed Fisher information
inv.fisher <- solve(nout$hessian)
se.inv.fisher <- sqrt(diag(inv.fisher))
se.inv.fisher

# to get standard error for gamma.hat we need to know
# var(a + b X) = b^2 Var(X) for constants a and b
# var(X + Y) = var(X) + var(Y) + 2 cov(X, Y)
# var(X - Y) = var(X) + var(Y) - 2 cov(X, Y)
# so var(1 - alpha.hat - beta.hat) = var(alpha.hat + beta.hat)
#     = var(alpha.hat) + var(beta.hat) + 2 cov(alpha.hat, beta.hat)
se.gamma.hat <- sqrt(inv.fisher[1,1] + inv.fisher[2,2] + 2 * inv.fisher[1,2])

foo <- cbind(c(alpha.hat, beta.hat, gamma.hat), c(se.inv.fisher, se.gamma.hat))
rownames(foo) <- c("alpha", "beta", "gamma")
colnames(foo) <- c("estimate", "std. err.")
round(foo, 3)

# and, of course, point estimate +/- 1.96 standard errors is an asymptotic
# 95% confidence interval
bar <- cbind(foo[,1] - qnorm(0.975) * foo[,2],
    foo[,1] + qnorm(0.975) * foo[,2])
rownames(bar) <- rownames(foo)
colnames(bar) <- c("lower", "upper")
round(bar, 3)