R devel - Feb 2000 - constrained optimiser doesn't obey constraints (PR#411)

---------- Forwarded message ----------
Date: Fri, 4 Feb 2000 01:58:44 +0000 (GMT)
From: Thomas Lumley <thomas@biostat.washington.edu>
To: r-bugs@biostat.ku.dk
Subject: constrained optimiser doesn't obey constraints


 The constrained optimiser is calling the function with parameter values
outside the constrained range

 In this example the box is from c(0,0,0,0) to c(1,1,.5,1). The function
prints out  the parameter vector when it is out of range. This gives:

[1] "x= -0.001" "x= 1"      "x= 0"      "x=
1"
[1] "x= 0"     "x= 1.001" "x= 0"     "x=
1"
[1] "x= 0"      "x= 1"      "x= -0.001" "x=
1"
[1] "x= 0"     "x= 1"     "x= 0"     "x=
1.001"
[1] "x= NaN" "x= NaN" "x= NaN" "x= NaN"
[1] "x= NaN" "x= NaN" "x= NaN" "x= NaN"
 repeated a number of times.

 For these starting value it eventually reports convergence, for
 some others it goes forever.

---Example code---

revise.mix.binbb_function(x,loss,no.inform){
    if (any(is.nan(x)) || any(x<0) || any(x>c(1,1,.5,1))){
        print(paste("x=",x))
        return(Inf)
    }
    eta_x[1]
    pi1_x[2]
    omega1_x[3]
    p0_x[4]
    log.term_0
    for (i in 1:length(no.inform)){
        if ( (loss[i]-1) >= 0) {
              first.prod_1
              for (r in 0:(loss[i]-1)){
                first.prod_first.prod*(pi1+r*omega1)
             };} else {first.prod_1;}
            if ( (no.inform[i]-loss[i]-1) >=0) {
        second.prod_1
        for (r in 0:(no.inform[i]-loss[i]-1)) {
          second.prod_second.prod*(1-pi1+r*omega1)
        };} else {second.prod_1;}

        third.prod_1
        for (r in 0:(no.inform[i]-1)){
          third.prod_third.prod*(1+r*omega1)
        }

        common.term<-choose(no.inform[i],loss[i])
        first.part_eta*common.term*first.prod*second.prod/third.prod

whole.term_eta*common.term*first.prod*second.prod/third.prod+(1-eta)*common.term*p0^loss[i]*(1-p0)^(no.inform[i]-loss[i])
        
        log.term_log.term+log(whole.term)
   
    }
    ## log.term is the log-likelihood.
    ## For nlminb we need to minimize -log-likelihood which is
    ## equivalent to maximizing the log-likelihood
    return(-log.term)
}


init<-c(0.001, 0.001, 0.001, 0.220)


"sim.ns" <-
c(17, 17, 15, 15, 16, 18, 18, 19, 17, 12, 18, 19, 19, 15, 20, 
17, 7, 15, 18, 15, 19, 19, 19, 21, 16, 17, 18, 19, 21, 12, 16, 
16, 21, 14, 16, 8, 18, 15)
"loss.dat" <-
structure(c(17, 13, 8, 9, 10, 4, 1, 7, 3, 4, 2, 6, 3, 4, 7, 3, 
1, 2, 5, 4, 3, 4, 2, 1, 1, 7, 4, 5, 3, 0, 4, 6, 4, 1, 6, 3, 2, 
3), .Names = c("V1", "V2", "V3", "V4",
"V5", "V6", "V7", "V8",
"V9", "V10", "V11", "V12",
"V13", "V14", "V15", "V16",
"V17",
"V18", "V19", "V20", "V21",
"V22", "V23", "V24", "V25",
"V26",
"V27", "V28", "V29", "V30",
"V31", "V32", "V33", "V34",
"V35",
"V36", "V37", "V38"))
 

optim(init,fn=revise.mix.binbb,
               lower=c(0.00,0.00,0.00,0.00),upper=c(1,1,.5,1),
            no.inform=sim.ns,loss=loss.dat,method="L-BFGS-B")

--please do not edit the information below--

Version:
 platform = i686-unknown-linux
 arch = i686
 os = linux
 system = i686, linux
 status = Under development (unstable)
 major = 0
 minor = 99.0
 year = 2000
 month = January
 day = 29
 language = R

Search Path:
 .GlobalEnv, Autoloads, package:base






-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

On Fri, 4 Feb 2000 thomas@biostat.washington.edu wrote:
> 
> 
> ---------- Forwarded message ----------
> Date: Fri, 4 Feb 2000 01:58:44 +0000 (GMT)
> From: Thomas Lumley <thomas@biostat.washington.edu>
> To: r-bugs@biostat.ku.dk
> Subject: constrained optimiser doesn't obey constraints
> 
> 
>  The constrained optimiser is calling the function with parameter values
> outside the constrained range

It seems to be happening in the numerical gradient computation. If the
current parameter value is on the boundary then one or other of the points
used for the finite difference will be over the boundary.

This explains why the values are over the boundary by a consistent small
amount.

	-thomas

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Your example seems not valid R (it gives me syntax errors, and with the use
of _ with no space around is completely unreadable to me: I suspect some
lines got split in the wrong place), so I can't run it.  However, I think
the problem is that you

1) Are not supplying derivatives
2) Have a non-differentiable function over the space.
3) Started too near the boundary (the finite-difference step is 1e-3,
by default).

and so the finite-difference code is not obeying the constraints.
Really you need either to supply derivatives or have a differentiable
function over a region a little larger than the constraints.

It should be possible to write the finite-difference code only to 
step into the region: do you want to do it?


-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._