[R] Optimize function in R: unable to find maximum of logistic function

[William Dunlap]

Note that f(x) returns exactly zero for all x above 1, hence its estimated
derivative is 0 everywhere in that region.   You are asking optimize to find
the non-zero part of a function which is 0 in more than 80% of its domain.

You can put a trace on f() to see where optimize() looks:
   > trace(f, quote(cat("k=", deparse(k), "\n")), print=FALSE)
  [1] "f"
  > optimize(f, c(0, 5), tol = 1e-14, maximum= TRUE)
  k= 1.90983005625053 
  k= 3.09016994374947 
  k= 3.81966011250105 
  k= 4.27050983124842 
  k= 4.54915028125263 
  k= 4.72135954999579 
  k= 4.82779073125683 
  k= 4.89356881873896
  ... eventually homing in on 5 ...
It looks in a few places, but f() is 0 in all of them so it figures that
is the mininum and the maximum value it ever takes.  I suppose
you could ask that optimize look in more places for a place with
a non-zero derivative but it would be very expensive to look in all
(c. 2^50) places.

You probably should have it maximize the log of your objective
function, which gives you more range to play with, and you will
have to do some numerical analysis to minimize underflow and
roundoff error.  E.g., one could use
    f1 <- function (k)  {
        # log(f(k))
        T_s <- 20
        log(2) - 2 * T_s * k - log1p(exp(-2 * T_s * k))
   }
where log1p(x) calculates log(1+x) more accurately than a computer's
'log' and '+' can.  (Remember that 1+10^-17 exactly equals 1 here.)

One can do better yet by using the logspace_add(logx, logy) function
available in the C API to R, but there is currently not an R-language interface
to that (or logspace_subtract(logx,logy)).  They calculate log(exp(logx) +- exp(logy))
with minimal roundoff error.

You could probably also use the built-in plogis() function, with its log=TRUE and
perhaps lower.tail=FALSE arguments.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of Rolf Turner
> Sent: Monday, October 28, 2013 2:01 PM
> To: Shantanu MULLICK
> Cc: r-help at r-project.org
> Subject: Re: [R] Optimize function in R: unable to find maximum of logistic function
> 
> 
> This could be described as a bug, perhaps.  Or it could be described as
> an indication that numerical optimization is inevitably tricky. Notice that
> if you narrow down your search interval from [0,5] to [0,0.5] you get the
> right answer:
> 
>  > optimize(f, c(0, 0.5), maximum= TRUE,tol=1e-10)
> $maximum
> [1] 4.192436e-11
> 
> $objective
> [1] 1
> 
> I guess there's a problem with finding a gradient that is effectively
> (numerically) zero when "k" is equal to 5.
> 
> 
>      cheers,
> 
>      Rolf Turner
> 
> 
> On 10/29/13 06:00, Shantanu MULLICK wrote:
> > Hello Everyone,
> >
> > I want to perform a 1-D optimization by using the optimize() function. I
> > want to find the maximum value of a "logistic" function. The optimize()
> > function gives the wrong result.
> >
> > My code:
> > f= function (k) {
> > T_s = 20
> > result = (2- 2/(1+ exp(-2*T_s*k)))
> > return(result)
> > }
> > optimize(f, c(0, 5), tol = 0.0000000000001, maximum= TRUE)
> >
> > The maximum value for the function happens at k=0, and the maximum value is
> > 1. Yet  the optimise function, says that the maximum value happens at k=
> > 4.9995, and the maximum value is 0.
> >
> 
> ______________________________________________
> 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.