Dear S+ users,
has anyone experienced problems with uniroot(), in Splus 4.5?
In my case, it locates a (non-existing) solution outside the
given search interval. The S-news database on Statlib does not
show anything of this sort.
In the sample problem below, a root is searched in the
inteval ( 0.02418690, 0.07753657), with an existing solution
near 0.05408. However uniroot comes up with the the answer
-0.6160564 (!)
The code below should reproduce the problem.
Possibly, the problem is due to the recursive call to uniroot().
However, R 0.64 produces more sensible answers.
Regards
Adelchi Azzalini
#---------------------------------------------------
f.dev <- function(mu,y,crit) f.deviance(mu,y,F)-crit
f.deviance <- function(mu, y, display=F){
g <- function(x,d) mean(d/(1+x*d))
epsilon<-(1e-8)
W<-rep(NA,length=length(mu))
for(i in 1:length(mu)) {
d <- (y-mu[i])
interv<- c(-1/max(d)+epsilon,-1/min(d)-epsilon)
root <- uniroot(g,interval=interv,d=d)
lambda<-root$root
W[i] <- 2*sum(log(1+lambda*d))
if(is.na(W[i])| W[i]<0) warning("W<0 or W==NA")
}
if(length(mu)>1 & display) plot(mu,W,type="l")
return(W)
}
#-----------------
# sample problem:
#
y <- c(0.8931869, 0.6805792, 0.8349237, rep(0,97))
interv <- c( 0.02418690, 0.07753657)
crit <- 2.705543
print(f.dev(interv,y,crit))
# [1] -2.705491 3.881472
print(f.dev(0.05408, y,crit))
# [1] -0.0001011083
uniroot(f.dev, interval=interv, y=y, crit=crit)$root
# [1] -0.6160564
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|