[R] Hessian in box-constraint problem - concern OPTIM function

[Cleber Nogueira Borges]

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