[R] Re: conf.int. for ecdfs {was "Two questions"}

[Martin Maechler]

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