From: Cleber Nogueira Borges <klebyn_at_yahoo.com.br>

Date: Tue, 24 Jun 2008 22:04:36 -0300

)

)

)

plot( xi, yi, pch=19, main=paste("Passo",i,"de",steps,sep=" "), cex=2 )

Sys.sleep(.11)

###file = paste("ISO_",(i+100),".png", sep="") ###savePlot(filename=file, type ="png", device=dev.cur() ) }

> Cleber Nogueira Borges wrote:

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 25 Jun 2008 - 01:11:37 GMT

Date: Tue, 24 Jun 2008 22:04:36 -0300

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 ______________________________________________R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 25 Jun 2008 - 01:11:37 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Wed 25 Jun 2008 - 19:31:50 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*