# [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,

>      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,
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)

## 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)
title(titulo)
}

> '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

```