[R] nlm, hessian, and derivatives in obj function?

[Thomas W Blackwell]

Jeff  -

The function  obj()  which you define below is just
a bit peculiar, since inside the function it assigns
attributes to an object 'obj' with the same name as
the *function* but which has not previously been
defined inside the function.  Is this really what
you intended ?  I'm not enough of a guru to figure
out what R will do with this syntax.

-  tom blackwell  -  u michigan medical school  -  ann arbor  -

On Fri, 17 Oct 2003, Jeff D. Hamann wrote:

> I've been working on a new package and I have a few questions regarding the
> behaviour of the nlm function. I've been (for better or worse) using the nlm
> function to fit a linear model without suppling the hessian or gradient
> attributes in the objective function. I'm curious as to why the nlm requires
> 31 iterations (for the linear model), and then it doesn't work when I try to
> add the derivative information. I know using nlm for a linear model isn't
> the "optimal" method, but I would like to make sure the parameter estimates
> and the se's are matching before I attempt more difficult problems.
>
> rm(list=ls(all=TRUE))
> print( "running nlsystemfit models test at end...")
> data( kmenta )
> attach( kmenta )
> ##demand2 <- q ~ d0 + d1 * p + d2 * d
> supply2 <- q ~ s0 + s1 * p + s2 * f + s3 * a
> ##system2 <- list( demand2, supply2 )
> ##labels <- list( "Demand", "Supply" )
> ##inst <- ~ d + f + a
> ##sv2 <- c(d0=3,s2=2.123,d2=4,s0=-2.123,s3=4.234,d1=4.234,s1=0.234)
> sv2 <- c(s0=-2.123,s1=0.234,s2=2.123,s3=4.234)
>
> obj <- function( s, eqn, data, parmnames )
> {
>   ## get the values of the parameters
>   for( i in 1:length( parmnames ) )
>     {
>       name <- names( parmnames )[i]
>       val <- s[i]
>       storage.mode( val ) <-  "double"
>       assign( name, val )
>     }
>
>   lhs <- as.matrix( eval( as.formula( eqn )[[2]] ) )
>   rhs <- as.matrix( eval( as.formula( eqn )[[3]] ) )
>   resid <- crossprod( lhs - rhs )
>
>   ## just how does this work...
>   attr( obj, "value" ) <- resid
>   attr( obj, "gradient" ) <- attr( eval( deriv3( eqn, names(
> parmnames ) ) ), "gradient" )
>
> }
>
> res <- nlm( obj, sv2, hessian=T, eqn=supply2, data=kmenta, parmnames=sv2,
> check.analyticals=T)
>
> I haven't been able to get nlm to function as I keep getting the following
> error message:
>
> Error in nlm(obj, sv2, hessian = T, eqn = supply2, data = kmenta, parmnames
> = sv2,  :
>  invalid function value in 'nlm' optimizer
>
>
>
> I was under the impression that you could also obtain the se of the
> parameter estimates using the sqrt( diag( res$hessian ) ), but I haven't
> been able to reproduce the se computed by the Jacobian
>
> se <- sqrt( mse * diag( solve( crossprod( J ) ) ) )    # gives the correct
> results...
> hse <- sqrt( ( res$minimum / 8 ) * diag( solve( res$hessian ) ) )   # gives
> similar results, but why 8?
>
> I've tried to put the functionality to include the jacobian and hessian in
> the objective function for nlm without success as I don't know what the form
> of the functions will be ahead of time.
>
> and get the se from the sqrt( diag( hessian ) ), but it's nowhere close?
>
> Jeff.
>
> ---
> Jeff D. Hamann
> Hamann, Donald and Associates, Inc.
> PO Box 1421
> Corvallis, Oregon USA 97339-1421
> (office) 541-754-1428
> (cell) 541-740-5988
> jeff_hamann at hamanndonald.com
> www.hamanndonald.com