[Rd] constrained optimiser doesn't obey constraints (PR#411)

[thomas@biostat.washington.edu]

---------- 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._