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 > [EMAIL PROTECTED] > www.hamanndonald.com ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help