[S] uniroot -- doesn't work recursively

[Martin Maechler]

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