[R] Hessian in box-constraint problem - concern OPTIM function
Cleber Nogueira Borges
klebyn at yahoo.com.br
Wed Jun 25 03:04:36 CEST 2008
Hello Mr. Graves,
Hello all useRs,
Many thanks for your attention.
> This is an interesting question.
> What is the problem you are trying to solve and how do the
> boundary conditions function as part of this system?
One of the functions that I need to minimized is:
sum( ( ( hat_xi - xi )^2 )*uxi^-2 + ( ( ( a+b*hat_xi ) - yi)^2 )*uyi^-2)
the context: I need to considerate errors in regressors: x[i] ~ N( x[i] ;
ux[i]^2 )
u = 'uncertainty' is the same of std, but this 'u' is because the
metrology terminology.
And, I would like to restrict the x[i] variables in ~95% CI.
My dirty code (test) follow below
Thanks again.
Cleber
###############
n <- 7 ### TODO: any number???
xi <- c(1:n) ; uxi <- round( abs( rnorm( n,0,1e-1)),6)
yi <- round(xi + uxi + rnorm(n,0,.9),6) ; uyi <- round(abs(rnorm(
n,0,1e-1)),6)
naive <- lm( yi ~ xi )
# p: parameters
p <- 2
plot( xi,yi )
abline( naive )
fobjetiva <- function( optPar=c(xi,a,b) , xi, uxi, yi, uyi )
{
n <- length( xi )
hat_xi <- optPar[1:n] ; a <- optPar[n+1] ; b <- optPar[n+2]
sum( ( (hat_xi - xi)^2 )*uxi^-2 + ( ( (a+b*hat_xi) - yi)^2
)*uyi^-2 )
}
## testing
fobjetiva(c(xi, coef(naive)), xi,uxi,yi,uyi)
### box-constraints
seCoef <- sqrt(diag(vcov( naive )))
Linf <- as.numeric(c( xi-2*uxi, coef(naive)-5000*seCoef ))
Lsup <- as.numeric(c( xi+2*uxi, coef(naive)+5000*seCoef ))
metodo <- 4 # 1,2,3,4 ou 5
all_iterOptim <- capture.output(
reportOptim <- optim(
par=c(xi,coef(naive)),
fn=fobjetiva,
xi=xi,
uxi=uxi,
yi=yi,
uyi=uyi,
method=c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN")[metodo],
lower=Linf,
upper=Lsup,
hessian=TRUE,
control=list(
REPORT=1,
maxit=1e4,
trace=6
)
)
)
### get all steps
all_steps <- grep("^X", all_iterOptim )
all_steps <- all_iterOptim[ all_steps ]
steps <- length(all_steps)
steps
matrix_steps <- matrix(0,nr=steps, nc=(n+p) )
for( i in 1:steps )
{
matrix_steps[i,] <- as.numeric(unlist(strsplit(all_steps[i], "
"))[3:(2+n+p)])
}
windows(restore=T)
### view a animation of this otimization
for( i in 1:steps )
{
x <- matrix_steps[i,1:n]
a <- matrix_steps[i,n+p-1]
b <- matrix_steps[i,n+p]
y <- a+b*x
plot( xi, yi, pch=19, main=paste("Passo",i,"de",steps,sep=" "), cex=2 )
abline(a,b, col='blue', lwd=3)
abline(naive, lwd=2, lty=2, col='red')
points( x, y, col='red', pch=19, cex=2 )
segments(xi,yi,x,y, lwd=2)
Sys.sleep(.11)
###file = paste("ISO_",(i+100),".png", sep="")
###savePlot(filename=file, type ="png", device=dev.cur() )
}
###comando = "convert -dispose previous -adjoin -delay 35 ISO_*.png
-loop 0 ISO_animator.gif"
###shell(comando)
###unlink("ISO_*.png")
## view trajectory
windows(restore=T)
par( mfrow=c(4,2))
for( i in 1:7 ){
restri <- xi[i]+uxi[i]*c(-2,2)
interv <- range(matrix_steps[,i], restri )
plot(1:steps, matrix_steps[,i], t='l', xlab="Passos", ylab="x1",
ylim=interv, lwd=2, las=2 )
abline( h=restri, col='red', lwd=2)
titulo=c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN")[metodo]
title(titulo)
}
> I ask because the asymptotic theory behind your formula for
> 's.e.' breaks down with parameters at boundaries. It assumes that you
> are minimizing the negative log(likelihood) AND the optimum is in the
> interior of the region AND the log(likelihood) is sufficiently close
> to being parabolic that a reasonable approximation for the
> distribution of the maximum likelihood estimates (MLEs) has a density
> adequately approximated by a second-order Taylor series expansion
> about the MLEs. In this case, transforming the parameters will not
> solve the problem. If the maximum is at a boundary and if you send
> the boundary to Inf with a transformation, then a second-order Taylor
> series expansion of the log(likelihood) about the MLEs will be locally
> flat in some direction(s), so the hessian can not be inverted.
> These days, the experts typically approach problems like this
> using Monte Carlo, often in the form of Markov Chain Monte Carlo
> (MCMC). One example of an analysis of this type of problem appears in
> section 2.4 of Pinheiro and Bates (2000) Mixed-Effects Models in S and
> S-Plus (Springer).
>
> Hope this helps. Spencer Graves
>
> Cleber Nogueira Borges wrote:
>>
>> Hello all useRs,
>>
>> I am using the OPTIM function with particular interest in the method
>> L-BFGS-B,
>> because it is a box-constraint method.
>> I have interest in the errors estimates too.
>> I make:
>> s.e. <- sqrt( diag( solve( optim(...,method='L-BFGS-B',
>> hessian=TRUE)$hessian )))
>> but in help say:
>> "Note that this is the Hessian of the unconstrained problem even if the
>> box constraints are active."
>> My doubts is:
>> How to obtain a authentic hessian for a box-constraint problem?
>> How I should make a interpretation of this result (concern the
>> hessian) ?
>> Is possible make some transformation or so can I considerate this
>> result a good approximation??
>> I am grateful for some help!
>> References are welcome! :-D
>> Cleber Borges
More information about the R-help
mailing list