R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

>  #  chex94rp.r            October 2007, 2008
>  # 
>  # ML analysis of Extended Example -- reparameterized
>  #
>  # first define likelihood function
>     ulglik <- function(alpha) { # log-likelihood for extended example
+                                 #  no gradient, Hessian
+         t1 <- exp(alpha[1])/(1+exp(alpha[1])+exp(alpha[2]))
+         t2 <- exp(alpha[2])/(1+exp(alpha[1])+exp(alpha[2]))
+         n <- c(17,182,60,176)
+         p <- c( 2*t1*t2, t1*(2-t1-2*t2), t2*(2-t2-2*t1), (1-t1-t2)**2 )
+         ulglik <- -crossprod( n, log(p) )   }
>  #  
>  #  use nlm for optimization
>   th1 <- nlm(ulglik, c(0,0), hessian=T, print.level=2)
iteration = 0
Step:
[1] 0 0
Parameter:
[1] 0 0
Function Value
[1] 678.145
Gradient:
[1]  30.33341 193.00009

iteration = 1
Step:
[1] -0.4134901 -2.6308825
Parameter:
[1] -0.4134901 -2.6308825
Function Value
[1] 529.771
Gradient:
[1]  87.71047 -42.93979

iteration = 2
Step:
[1] -0.963755  0.090721
Parameter:
[1] -1.377245 -2.540161
Function Value
[1] 515.0408
Gradient:
[1] -54.50460 -27.73849

iteration = 3
Step:
[1] 0.3954172 0.4092844
Parameter:
[1] -0.981828 -2.130877
Function Value
[1] 494.122
Gradient:
[1]  -9.463442 -11.190854

iteration = 4
Step:
[1] 0.1056419 0.2153792
Parameter:
[1] -0.8761862 -1.9154979
Function Value
[1] 492.549
Gradient:
[1] 1.3045229 0.8245283

iteration = 5
Step:
[1] -0.01136476 -0.01606467
Parameter:
[1] -0.887551 -1.931563
Function Value
[1] 492.5353
Gradient:
[1]  0.01819541 -0.05483036

iteration = 6
Step:
[1] -1.903352e-05  7.943617e-04
Parameter:
[1] -0.887570 -1.930768
Function Value
[1] 492.5353
Gradient:
[1] -0.0015572255  0.0009447857

iteration = 7
Parameter:
[1] -0.8875606 -1.9307788
Function Value
[1] 492.5353
Gradient:
[1] 1.358558e-05 8.832200e-08

Relative gradient close to zero.
Current iterate is probably solution.

>   #
>   # now get inverse of Hessian for std errors
>       th1$hessian
          [,1]      [,2]
[1,] 143.46939 -21.43636
[2,] -21.43636  69.72793
>       L <- t( chol(th1$hessian) )
>       L
          [,1]     [,2]
[1,] 11.977871 0.000000
[2,] -1.789664 8.156288
>       backsolve( t(L), forwardsolve(L,diag(1,nrow(L))) ) 
            [,1]        [,2]
[1,] 0.007305710 0.002245984
[2,] 0.002245984 0.015031934
>   # 
>   # probabilities in the original parameterization
>      exp(th1$estimate[1])/(1+exp(th1$estimate[1])+exp(th1$estimate[2]))
[1] 0.2644442
>      exp(th1$estimate[2])/(1+exp(th1$estimate[1])+exp(th1$estimate[2]))
[1] 0.09316873
>   #  done
>     rm(list=ls())
>   q()