[R] Re: conf.int. for ecdfs {was "Two questions"}
Martin Maechler
maechler at stat.math.ethz.ch
Mon Oct 22 19:15:35 CEST 2001
>>>>> "Christian" == Christian Blouin <bongo at havana.biochem.dal.ca> writes:
Christian> I have two questions that I could not answer from the
Christian> documentation.
Christian> A - ecdf and confidence intervals : Is there a (simple) way
Christian> to generate confidence intervals (95%) for a ecdf?
About two weeks ago, Kjetil Halvorsen (CC'ed above)
sent me code to just do this, and I have improved on the code a bit
with the idea to incorporate into the stepfun package (part of "standard R").
However, it's content is also linked with the (hidden) pkstwo() function
in ctest's kolomogorov-smirnov testing.
Also, I think part of this might become an option to ecdf() or plot.ecdf()
rather than the current setup.
Unfortunately, it's not finished for integration into R,
but here is the current code (without help pages) that works quite nicely:
## Fixme: In the following, computing and plotting should be separated
###--> ./ecdf.R plot.ecdf() should get conf.type and conf.int argument!!
ecdf.ksCI <- function(x, main = NULL, sub = NULL,
xlab = deparse(substitute(x)), ...)
{
xlab
if(is.null(main))
main <- paste("ecdf(",deparse(substitute(x)),") + 95% K.S. bands",
sep="")
n <- length(x)
if(is.null(sub))
sub <- paste("n = ", n)
ec <- ecdf(x)
xx <- get("x", envir=environment(ec))# = sort(x)
yy <- get("y", envir=environment(ec))
D <- approx.ksD(n)
yyu <- pmin(yy+D, 1)
yyl <- pmax(yy-D, 0)
ecu <- stepfun(xx, c(yyu, 1) )
ecl <- stepfun(xx, c(yyl, yyl[n]) )
## Plots -- all calling plot.stepfun
plot(ec, main = main, sub = sub, xlab = xlab)## no "..." here ??? __hmm..__
plot(ecu, add=TRUE, verticals=TRUE, do.points=FALSE,
col.hor="red" , col.vert="red", ...)
plot(ecl, add=TRUE, verticals=TRUE, do.points=FALSE,
col.hor="red", col.vert="red", ...)
}
approx.ksD <- function(n)
{
## approximations for the critical level for Kolmogorov-Smirnov statistic D,
## for confidence level 0.95. Taken from Bickel & Doksum, table IX, p.483
## and Lienert G.A.(1975) who attributes to Miller,L.H.(1956), JASA
ifelse(n > 80,
1.358 /( sqrt(n) + .12 + .11/sqrt(n)),##Bickel&Doksum, table IX,p.483
splinefun(c(1:9, 10, 15, 10 * 2:8),# from Lienert
c(.975, .84189, .70760, .62394, .56328,# 1:5
.51926, .48342, .45427, .43001, .40925,# 6:10
.33760, .29408, .24170, .21012,# 15,20,30,40
.18841, .17231, .15975, .14960)) (n))
}
###-- Example from Kjetil :
ecdf.ksCI( rchisq(50,3) )
--
Martin Maechler <maechler at stat.math.ethz.ch> http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum LEO D10 Leonhardstr. 27
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-1-632-3408 fax: ...-1228 <><
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help 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-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list