#-*- R -*-## Script from Fourth Edition of `Modern Applied Statistics with S'# Chapter 16   Optimization and Mazimum Likelihood Estimationlibrary(MASS)postscript(file="ch16.ps", width=8, height=8, pointsize=9)options(echo = T, width=65, digits=5)# 16.3 General optimizationattach(geyser)truehist(waiting, xlim = c(35, 110), ymax = 0.04, h = 5)wait.dns <- density(waiting, n = 512, width = "SJ")lines(wait.dns, lty = 2)lmix2 <- deriv3(     ~ -log(p*dnorm((x-u1)/s1)/s1 + (1-p)*dnorm((x-u2)/s2)/s2),     c("p", "u1", "s1", "u2", "s2"),     function(x, p, u1, s1, u2, s2) NULL)(p0 <- c(p = mean(waiting < 70), u1 = 50, s1 = 5, u2 = 80, s2 = 5))## using optimmix.obj <- function(p, x){  e <- p[1] * dnorm((x - p[2])/p[3])/p[3] +       (1 - p[1]) * dnorm((x - p[4])/p[5])/p[5]  if(any(e <= 0)) Inf else -sum(log(e))}optim(p0, mix.obj, x = waiting)$par # Nelder-Meadoptim(p0, mix.obj, x = waiting, method = "BFGS",     control = list(parscale= c(0.1, rep(1, 4))))$par# with derivativeslmix2a <- deriv(     ~ -log(p*dnorm((x-u1)/s1)/s1 + (1-p)*dnorm((x-u2)/s2)/s2),     c("p", "u1", "s1", "u2", "s2"),     function(x, p, u1, s1, u2, s2) NULL)mix.gr <- function(p, x) {   u1 <- p[2]; s1 <- p[3]; u2 <- p[4]; s2 <- p[5]; p <- p[1]   colSums(attr(lmix2a(x, p, u1, s1, u2, s2), "gradient")) }optim(p0, mix.obj, mix.gr, x = waiting, method = "BFGS",     control = list(parscale= c(0.1, rep(1, 4))))$parmix.nl0 <- optim(p0, mix.obj, mix.gr, method = "L-BFGS-B", hessian = T,                lower = c(0, -Inf, 0, -Inf, 0),                upper = c(1, rep(Inf, 4)), x = waiting)rbind(est = mix.nl0$par, se = sqrt(diag(solve(mix.nl0$hessian))))dmix2 <- function(x, p, u1, s1, u2, s2)             p * dnorm(x, u1, s1) + (1-p) * dnorm(x, u2, s2)attach(as.list(mix.nl0$par))wait.fdns <- list(x = wait.dns$x,                  y = dmix2(wait.dns$x, p, u1, s1, u2, s2))lines(wait.fdns)par(usr = c(0, 1, 0, 1))legend(0.1, 0.9, c("Normal mixture", "Nonparametric"),       lty = c(1, 2), bty = "n")pmix2 <- deriv(~ p*pnorm((x-u1)/s1) + (1-p)*pnorm((x-u2)/s2),               "x", function(x, p, u1, s1, u2, s2) {})pr0 <- (seq(along = waiting) - 0.5)/length(waiting)x0 <- x1 <- as.vector(sort(waiting)) ; del <- 1; i <- 0while((i <- 1 + 1) < 10 && abs(del) > 0.0005) {  pr <- pmix2(x0, p, u1, s1, u2, s2)  del <- (pr - pr0)/attr(pr, "gradient")  x0 <- x0 - 0.5*del  cat(format(del <- max(abs(del))), "\n")}detach()par(pty = "s")plot(x0, x1, xlim = range(x0, x1), ylim = range(x0, x1),     xlab = "Model quantiles", ylab = "Waiting time")abline(0, 1)par(pty = "m")mix1.obj <- function(p, x, y){  q <- exp(p[1] + p[2]*y)  q <- q/(1 + q)  e <- q * dnorm((x - p[3])/p[4])/p[4] +       (1 - q) * dnorm((x - p[5])/p[6])/p[6]  if(any(e <= 0)) Inf else -sum(log(e))}p1 <- mix.nl0$par; tmp <- as.vector(p1[1])p2 <- c(a = log(tmp/(1-tmp)), b = 0, p1[-1])mix.nl1 <- optim(p2, mix1.obj, method = "L-BFGS-B",                lower = c(-Inf, -Inf, -Inf, 0, -Inf, 0),                upper = rep(Inf, 6), hessian = T,                x = waiting[-1], y = duration[-299])rbind(est = mix.nl1$par, se = sqrt(diag(solve(mix.nl1$hessian))))if(!exists("bwt")) {  attach(birthwt)  race <- factor(race, labels=c("white", "black", "other"))  ptd <- factor(ptl > 0)  ftv <- factor(ftv); levels(ftv)[-(1:2)] <- "2+"  bwt <- data.frame(low=factor(low), age, lwt, race,	   smoke=(smoke>0), ptd, ht=(ht>0), ui=(ui>0), ftv)  detach(); rm(race, ptd, ftv)}logitreg <- function(x, y, wt = rep(1, length(y)),               intercept = T, start = rep(0, p), ...){  fmin <- function(beta, X, y, w) {      p <- plogis(X %*% beta)      -sum(2 * w * ifelse(y, log(p), log(1-p)))  }  gmin <- function(beta, X, y, w) {      eta <- X %*% beta; p <- plogis(eta)      -2 * matrix(w *dlogis(eta) * ifelse(y, 1/p, -1/(1-p)), 1) %*% X  }  if(is.null(dim(x))) dim(x) <- c(length(x), 1)  dn <- dimnames(x)[[2]]  if(!length(dn)) dn <- paste("Var", 1:ncol(x), sep="")  p <- ncol(x) + intercept  if(intercept) {x <- cbind(1, x); dn <- c("(Intercept)", dn)}  if(is.factor(y)) y <- (unclass(y) != 1)  fit <- optim(start, fmin, gmin, X = x, y = y, w = wt,               method = "BFGS", ...)  names(fit$par) <- dn  cat("\nCoefficients:\n"); print(fit$par)  # R: use fit$value and fit$convergence  cat("\nResidual Deviance:", format(fit$value), "\n")  if(fit$convergence > 0)      cat("\nConvergence code:", fit$convergence, "\n")  invisible(fit)}options(contrasts = c("contr.treatment", "contr.poly"))X <- model.matrix(terms(low ~ ., data=bwt), data = bwt)[, -1]logitreg(X, bwt$low)AIDSfit <- function(y, z, start=rep(mean(y), ncol(z)), ...){  deviance <- function(beta, y, z) {      mu <- z %*% beta      2 * sum(mu - y - y*log(mu/y)) }  grad <- function(beta, y, z) {      mu <- z %*% beta      2 * t(1 - y/mu) %*% z }  optim(start, deviance, grad, lower = 0, y = y, z = z,        method = "L-BFGS-B", ...)}Y <- scan()12 14 33 50 67 74 123 141 165 204 253 246 240library(nnet) # for class.inds <- seq(0, 13.999, 0.01); tint <- 1:14X <- expand.grid(s, tint)Z <- matrix(pweibull(pmax(X[,2] - X[,1],0), 2.5, 10),length(s))Z <- Z[,2:14] - Z[,1:13]Z <- t(Z) %*% class.ind(factor(floor(s/2))) * 0.01round(AIDSfit(Y, Z)$par)rm(s, X, Y, Z)# End of ch16