R devel - Apr 1999 - [S] uniroot -- doesn't work recursively

Dear Prof Azzalini,

You have an interesting example of  recursive use of  uniroot().
    [re-cited at the end]
However, note that R currently has the same problem as S-plus:

Uniroot() doesn't work reliably, recursively.
When you found it to be better, you were just lucky.

The relevant file in R's source,   src/main/optimize.c
says 

   /* WARNING : As things stand, these routines should not be called
    *	     recursively because of the way global variables are used.
    *	     This could be fixed by saving and restoring these global variables.
    */

--------

Fellow R-code developpers :  

       Does anyone of you feel like fixing these?
       Would be really nice and useful!

       (and you will become famous in the R community ... :-)

Martin Maechler <maechler@stat.math.ethz.ch>
http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum SOL G1;	Sonneggstr.33
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1086			<><
>>>>> "AA" == Adelchi Azzalini
<aa@gwen.stat.unipd.it> writes
      to S-news on 22 Apr 1999 :

    AA> Dear S+ users,

    AA> has anyone experienced problems with uniroot(), in Splus 4.5?  
    AA> In my case, it locates a (non-existing) solution outside the 
    AA> given search interval. The S-news database on Statlib does not
    AA> show anything of this sort.

    AA> In the sample problem below, a root is searched in  the 
    AA> interval ( 0.02418690, 0.07753657), with an existing solution
    AA> near  0.05408. 
    AA> However uniroot comes up with the the answer -0.6160564 (!)

    AA> The code below should reproduce the problem.
    AA> Possibly, the problem is due to the recursive call to uniroot().
    AA> However, R 0.64 produces more sensible answers.

    AA> Regards
    AA> Adelchi Azzalini

    AA> #---------------------------------------------------

    AA> f.dev <- function(mu,y,crit) f.deviance(mu,y,F)-crit

    AA> f.deviance <- function(mu, y, display=F){
    AA> g <- function(x,d) mean(d/(1+x*d))
    AA> epsilon<-(1e-8)
    AA> W<-rep(NA,length=length(mu))
    AA> for(i in 1:length(mu)) {
    AA> d <- (y-mu[i])
    AA> interv<- c(-1/max(d)+epsilon,-1/min(d)-epsilon)
    AA> root <- uniroot(g,interval=interv,d=d) 
    AA> lambda<-root$root
    AA> W[i] <- 2*sum(log(1+lambda*d))
    AA> if(is.na(W[i])| W[i]<0) warning("W<0 or W==NA")
    AA> }
    AA> if(length(mu)>1 & display) plot(mu,W,type="l")
    AA> return(W)
    AA> }

    AA> #-----------------
    AA> #  sample problem:
    AA> #
    AA> y <- c(0.8931869, 0.6805792, 0.8349237, rep(0,97))
    AA> interv <- c( 0.02418690, 0.07753657)
    AA> crit <- 2.705543

    AA> print(f.dev(interv,y,crit))
    AA> # [1] -2.705491  3.881472

    AA> print(f.dev(0.05408, y,crit))
    AA> # [1] -0.0001011083

    AA> uniroot(f.dev, interval=interv, y=y, crit=crit)$root
    AA> # [1] -0.6160564
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Martin Maechler <maechler@stat.math.ethz.ch> writes:
> Dear Prof Azzalini,
> 
> You have an interesting example of  recursive use of  uniroot().
...
> Fellow R-code developpers :  
> 
>        Does anyone of you feel like fixing these?
>        Would be really nice and useful!
> 
>        (and you will become famous in the R community ... :-)
Well, uniroot calls zeroin which is a FORTRAN routine. As long as this
is the case, all bets are off as to whether or not it will be
reentrant on any given system. In principle, it is not, and FORTRAN
just plain states that "you cannot do that", but if your compiler
happens to be f2c-like, and no global storage is involved, then it
*might* work after all. The example doesn't work on my system either.

The only sensible thing to do to get this to work is to rewite zeroin
in C. According to the comments in the file, it is a transcription of
an Algol60 program, so it is probably not hard to convert it to
another Algol descendant.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._