EXAMINE THE EVALUATION OF MAXIMA FOR LIKELIHOOD-LIKE 
FUNCTIONS WITH AND WITHOUT GRADIENT AND HESSIAN ATTRIBUTES
=======================================================

### Simple case:  2500 Weibull with shape 1.6, scale 0.4

xdat = rweibull(2500,1.6,0.4)

logLWeib0 = function (par2,dvec) {
      n = length(dvec)
      a = par2[1]
      b = par2[2]
      -log(a)+a*log(b) - (a-1)*mean(log(dvec)) + 
           mean(dvec^a)/b^a }

unix.time({ tmp = nlm(logLWeib0, c(1,1), dvec=xdat) })
   user  system elapsed 
   0.05    0.00    0.04 
> str(tmp)
List of 5
 $ minimum   : num -0.145
 $ estimate  : num [1:2] 1.559 0.403
 $ gradient  : num [1:2] -2.39e-07 1.14e-07
 $ code      : int 1
 $ iterations: int 8
## gradients of order 1e-7: NB accurate parameter est's!

logLWeib1 = function (par2,dvec) {
      n = length(dvec)
      a = par2[1]
      b = par2[2]
      obj = -log(a)+a*log(b) - (a-1)*mean(log(dvec)) + 
           mean(dvec^a)/b^a
      attr(obj,"gradient") = c(-1/a + log(b) - 
            mean(log(dvec)) + mean(log(dvec/b)*dvec^a)/b^a,
            a/b- a*mean(dvec^a)/b^(a+1))   
      obj }

> logLWeib0(c(1,1),xdat)
[1]  0.3620798
> logLWeib1(c(1,1),xdat)
[1]  0.3620798
attr(,"gradient")
[1] -0.009937885  0.637920180

> (logLWeib0(c(1+1e-6,1),xdat)-logLWeib0(c(1-1e-6,1),xdat))/2e-6
[1] -0.009937885
> (logLWeib0(c(1,1+1e-6),xdat)-logLWeib0(c(1,1-1e-6),xdat))/2e-6
[1] 0.6379202

unix.time({ tmp = nlm(logLWeib1, c(1,1), dvec=xdat) })
   user  system elapsed 
   0.07    0.00    0.07 
### same result but longer computation ! 

unix.time({ tmp = nlm(logLWeib0, c(1.6,0.4), dvec=xdat) })
   user  system elapsed 
   0.06    0.00    0.06   ## Note longer time, more iterations !
> str(tmp)
List of 5
 $ minimum   : num -0.145
 $ estimate  : num [1:2] 1.559 0.403
 $ gradient  : num [1:2] -4.17e-08 3.61e-07
 $ code      : int 1
 $ iterations: int 10


logLWeib2 = function (par2,dvec) {
      n = length(dvec)
      a = par2[1]
      b = par2[2]
      obj = -log(a)+a*log(b) - (a-1)*mean(log(dvec)) + 
           mean(dvec^a)/b^a
         attr(obj,"gradient") = c(-1/a + log(b) - 
            mean(log(dvec)) + mean(log(dvec/b)*dvec^a)/b^a,
            a/b- a*mean(dvec^a)/b^(a+1))   
      attr(obj,"hessian") = array(c(1/a^2 + 
            mean(log(dvec/b)^2*dvec^a)/b^a, rep(1/b-
            mean(dvec^a)/b^(a+1)-(a/b^(a+1))*mean(log(dvec/b)*
            dvec^a),2), -a/b^2 + (a*(a+1)/b^(a+2))*mean(dvec^a)),
            c(2,2))
      obj }


> logLWeib2(c(1.6,0.4),xdat)
[1] -0.1445586
attr(,"gradient")
[1]  0.03327704 -0.09032663
attr(,"hessian")
          [,1]      [,2]
[1,]  0.732797 -1.228535
[2,] -1.228535 16.587123


> (attr(logLWeib1(c(1.6+1e-6,.4),xdat),"gradient")-
+  attr(logLWeib1(c(1.6-1e-6,.4),xdat),"gradient"))/2e-6
[1]  0.732797 -1.228535
> (attr(logLWeib1(c(1.6,.4+1e-6),xdat),"gradient")-
+  attr(logLWeib1(c(1.6,.4-1e-6),xdat),"gradient"))/2e-6
[1] -1.228535 16.587123


> logLWeib2(c(1,1),xdat)
[1] 0.3620798
attr(,"gradient")
[1] -0.009937885  0.637920180
attr(,"hessian")
         [,1]       [,2]
[1,] 1.357504  1.2914855
[2,] 1.291486 -0.2758404     ### Note not pos-def at (1,1) !

unix.time({ tmp = nlm(logLWeib2, c(1,1), hessian=T, dvec=xdat) })
   user  system elapsed 
   0.29    0.00    0.29 
> str(tmp)
List of 6
 $ minimum   : num -0.145
 $ estimate  : num [1:2] 1.559 0.403
 $ gradient  : num [1:2] 1.88e-12 -1.75e-07
 $ hessian   : num [1:2, 1:2] 0.741 -1.041 -1.041 14.966
 $ code      : int 1
 $ iterations: int 8
## Note more accurate answer , gradient smaller

attr(logLWeib2(c(tmp$est),xdat),"hessian")
          [,1]      [,2]
[1,]  0.740968 -1.041975
[2,] -1.041975 14.983215
tmp$hessian
           [,1]      [,2]
[1,]  0.7409032 -1.040887
[2,] -1.0408871 14.966268    ## NOTE not exactly the same !

### One check on the validity of the model assumption 
### is to compare the Fisher Info found here with the
### sum of gradients time their transposes: the coding 
### is as follows:

> WeibScor = function(par2,dvec) {
      n = length(dvec)
      a = par2[1]
      b = par2[2]
      scor = cbind(-1/a + log(b) - log(dvec) + 
           log(dvec/b)*(dvec/b)^a, a/b- (a/b^(a+1))*dvec^a)
      list(scor=scor,  gradllik = c(t(scor) %*% rep(1,n)),
           altInfo = t(scor) %*% scor/n) }



> xdat = rweibull(2500,1.6,0.4)
> logLWeib2(c(1.6,0.4),xdat)
[1] -0.1710576
attr(,"gradient")
[1] -0.004228303 -0.005230492
attr(,"hessian")
          [,1]      [,2]
[1,]  0.713670 -1.060161
[2,] -1.060161 16.033998

> unlist(WeibScor(c(1.6,0.4),xdat)[2:3])
  gradllik1   gradllik2    altInfo1    altInfo2    altInfo3    altInfo4
-10.5707566 -13.0762305   0.7369663  -1.1290419  -1.1290419  16.0853598