thr3ads.net - R help - [R] nlm, hessian, and derivatives in obj function? [Oct 2003]

If this information is useful, please help other people find it:
Share via:

Jeff D. Hamann

2003-Oct-17 23:54 UTC

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

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


If I perform the fit without the derivative information, I get the correct
estimates,

$minimum
[1] 92.55106

$estimate
[1] 58.2754312  0.1603666  0.2481333  0.2483023

$gradient
[1] 8.552542e-08 9.087699e-06 5.716032e-06 2.163105e-06

$hessian
         [,1]       [,2]     [,3]     [,4]
[1,]   40.000   4000.762   3865.0   420.00
[2,] 4000.762 401486.918 386045.8 42007.76
[3,] 3865.000 386045.812 379593.1 39762.40
[4,]  420.000  42007.764  39762.4  5740.00

$code
[1] 1

$iterations
[1] 31

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

Thomas W Blackwell

2003-Oct-18 19:37 UTC

head link

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

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

Arne Henningsen

2003-Oct-23 11:05 UTC

head link

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

Hi,

I don't know much about non-linear models, but there is another possibility
to
fit these models: 

1) get some starting values for the parameters
2) take the derivatives of the model with respect to the parameters at the 
point of the starting values of the parameters
3) perform a linear estimation of this linearized model (using systemfit) to 
get new parameter estimates
4) got to step 2) and take these new parameter estimates in place of the 
starting values
5) iterate this until the parameters stay stable from one to the next 
iteration

This has three advantages:
1) It is not much work to write these function since systemfit already exists
2) If the model is linear in parameters, it is identical to the linearized 
model and, thus, the first iteration leads directly to the optimum
3) You get get the SEs from the last iteration of systemfit

Does this approach also have disadvantages (e.g. non-convergence of parameters 
in many cases)?

Best wishes,
Arne


On Saturday 18 October 2003 01:54, 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
>
>
> If I perform the fit without the derivative information, I get the correct
> estimates,
>
> $minimum
> [1] 92.55106
>
> $estimate
> [1] 58.2754312  0.1603666  0.2481333  0.2483023
>
> $gradient
> [1] 8.552542e-08 9.087699e-06 5.716032e-06 2.163105e-06
>
> $hessian
>          [,1]       [,2]     [,3]     [,4]
> [1,]   40.000   4000.762   3865.0   420.00
> [2,] 4000.762 401486.918 386045.8 42007.76
> [3,] 3865.000 386045.812 379593.1 39762.40
> [4,]  420.000  42007.764  39762.4  5740.00
>
> $code
> [1] 1
>
> $iterations
> [1] 31
>
> 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
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
-- 
Arne Henningsen
Department of Agricultural Economics
Christian-Albrechts-University Kiel
24098 Kiel, Germany
Tel: +49-431-880-4445
Fax: +49-431-880-1397 
ahenningsen at email.uni-kiel.de
http://www.uni-kiel.de/agrarpol/ahenningsen/

Possibly Parallel Threads

Search for more apparently analagous threads

R help - Oct 2003 - nlm, hessian, and derivatives in obj function?

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

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

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

Possibly Parallel Threads