Homework 13 solution, March 2008.
=================================

> HW13Data = rlogis(5000, runif(1,2,2.5), runif(1,.7,.85))
> mean(HW13Data)
[1] 2.480233
> sqrt(var(HW13Data))
[1] 1.366586
> integrate(function(x) x^2*dlogis(x), -20,20,abs.tol=1e-9)$val
[1] 3.289866

> unix.time({xfit = nlm(function(x) 
           -sum(log(dlogis(HW13Data,x[1],x[2]))), c(2,1))})
[1] 0.12 0.00 0.13   NA   NA
> str(xfit)
List of 5
 $ minimum   : num 8573
 $ estimate  : num [1:2] 2.48 0.75
 $ gradient  : num [1:2] -0.00251  0.00363
 $ code      : int 1
 $ iterations: int 7

> unix.time({xfit = nlm(function(x) 
           -sum(log(dlogis(HW13Data,x[1],x[2]))), c(0,1))})
[1] 0.19 0.00 0.20   NA   NA

> unix.time({xfit = nlm(function(x) -sum(log(dlogis(
           HW13Data,x[1],x[2]))), c(0,1), hessian=T)})
[1] 0.18 0.02 0.20   NA   NA

> unix.time({xfit2 = nlm(function(x) { 
         obj = -sum(log(dlogis(HW13Data,x[1],x[2])))
         pvec=plogis(HW13Data,x[1],x[2])
         nd = length(HW13Data)
         attr(obj,"gradient") = c(nd - 2*sum(pvec), 
              nd-sum((HW13Data-x[1])*(2*pvec-1))/x[2])/x[2]
              obj}, c(2,1))})
[1] 0.15 0.00 0.15   NA   NA
> str(xfit2)
List of 5
 $ minimum   : num 8573
 $ estimate  : num [1:2] 2.48 0.75
 $ gradient  : num [1:2] -0.00251  0.00364
 $ code      : int 1
 $ iterations: int 7
> unix.time({xfit2 = nlm(function(x) { 
         obj = -sum(log(dlogis(HW13Data,x[1],x[2])))
         pvec=plogis(HW13Data,x[1],x[2])
         nd = length(HW13Data)
         attr(obj,"gradient") = c(nd - 2*sum(pvec), 
              nd-sum((HW13Data-x[1])*(2*pvec-1))/x[2])/x[2]
              obj}, c(2,1), hessian=T)})
[1] 0.18 0.00 0.17   NA   NA

> unix.time({xfit2 = nlm(function(x) { 
         obj = -sum(log(dlogis(HW13Data,x[1],x[2])))
         pvec=plogis(HW13Data,x[1],x[2])
         nd = length(HW13Data)
         attr(obj,"gradient") = c(nd - 2*sum(pvec), 
              nd-sum((HW13Data-x[1])*(2*pvec-1))/x[2])/x[2]
              obj}, c(0,1))})
[1] 0.16 0.00 0.16 

> unix.time({xfit3 = nlm(function(x) { 
         dvec=dlogis(HW13Data,x[1],x[2])
         pvec=plogis(HW13Data,x[1],x[2])
         stdvec=(HW13Data-x[1])/x[2]
         obj = -sum(log(dvec))
         nd = length(HW13Data)
         attr(obj,"gradient") = c(nd - 2*sum(pvec), 
         nd-sum(stdvec*(2*pvec-1)))/x[2]
         attr(obj,"hessian") = -matrix(c(-2*sum(dvec),
           rep(nd-2*sum(pvec+stdvec*dvec),2), nd+2*sum(
           stdvec*(1-2*pvec-dvec*stdvec))), ncol=2)/x[2]^2
              obj}, c(2,1))})
[1] 0.74 0.00 0.73  NA NA

## So in simple problems or maybe not-so-simple there
##   might be no benefit in using the analytical 
##   gradients unless you want the output Hessian 
##   (which we almost always do want for statistical 
##   applications), but even then the analytical 
##   Hessian might be more trouble than it is worth !!