R/psistat.R

In GeoBosh/psistat: Statistical Functions for Course Unit Practical Statistics I at the University of Manchester

## Do not edit this file manually.
## It has been automatically generated from *.org sources.

psi.rbern <- function(n, p){
    rbinom(n, size = 1, prob = p)
}

# see ks.c in R sources
#  "pkolmogorov2x" computes the exact K_n(x) (cdf of Kolmogorov dist), e.g.
#      .C("pkolmogorov2x", p = as.double(x), n = as.integer(n), PACKAGE = "stats")$p
#
#  "pkstwo" computes the asymptotic approx L(d) such that K_n(x) approx=L(x*sqrt(n))
#         .C("pkstwo", as.integer(length(x)),
#                      p = as.double(x), as.double(tol), PACKAGE = "stats")$p
#
#  vizh sasto psmirnov2x()
#
# vizh kak se vikat tezi f. v ks.test.
# vizh i init.c v source-ovete.
#
# !!! Important difference: pkstwo accepts a vector x
#     while "pkolmogorov2x" computes K_n(x) for a single value of x only.

psi.pks <- function(q,n, lower.tail = TRUE, log.p=FALSE, exact=NULL,warn=FALSE){
    if(n<1)
        stop("n must be greater than or equal to 1.")

    if(is.null(exact)){
        if(n <= 100){
            exact <- TRUE
        }else{
            if(warn){
                print("n too large, using asymptotic approximation.")
                print("You may use argument exact=TRUE to force exact method for large n,")
                print("  but then the calculation may be very, very slow...")
            }
            exact <- FALSE
        }
    }

    if(n>100 && isTRUE(exact))
        print("Computing..., when n>100 the calculation may take forever.")

    f <- function(x){
        ## 2014-03-17
        ## .C("pkolmogorov2x", p = as.double(x), n = as.integer(n), PACKAGE = "stats")$p
        ##    ; also removing $p
        ## .Call(stats:::C_pKolmogorov2x,
        ##       p = as.double(x), n = as.integer(n), PACKAGE = "stats")
        ##
        ## copying the C functions in psistat, for CRAN
        check2 <- .Call(C_pKolmogorov2x,  p = as.double(x), n = as.integer(n))
                # # TODO: for testing:
                #     check1 <- .Call(stats:::C_pKolmogorov2x,
                #               p = as.double(x), n = as.integer(n), PACKAGE = "stats")
                #  if(check1 != check2)
                #      warning("results of stats:::C_pKolmogorov2x and C_pKolmogorov2x are not the same")
        check2
        }

    fasy <- function(x,tol=1e-06){
                # 2014-03-17
                #    ; also removing $p
                # .C("pkstwo", as.integer(length(x)),
                #   p = as.double(x), as.double(tol), PACKAGE = "stats")$p
                #    if(length(IND)) p[IND] <- .Call(C_pKS2, p = x[IND], tol)
                #
                ##   ks.c says: /* Two-sample two-sided asymptotic distribution */
                ##   but it seems to be the same for one-sample two-sided
                ##
                ## .Call(stats:::C_pKS2, # as.integer(length(x)),
                ##      p = x, tol)

                 # check1 <- .Call(stats:::C_pKS2, # as.integer(length(x)),
                 #              p = x, tol,  PACKAGE = "stats")

                check2 <- .Call(C_pKS2, # as.integer(length(x)),
               p = x, tol,  PACKAGE = "psistat")

                 # if(check1 != check2)
                 #     warning("results of stats:::C_pKS2 and C_pKS2 are not the same")
        check2
               ## todo: put stopifnot(check1 == check2)
           }


    wrk <- q
    ind <- 0<q & q < 1
    if(any(q<=0))
        wrk[q<=0] <- 0

    if(any(q>=1))
        wrk[q>=1] <- 1
    wrk[ind] <- if(isTRUE(exact)) sapply(q[ind],f) else fasy(sqrt(n)*q[ind])

    if(!lower.tail)
        wrk <- 1-wrk

    if(log.p){
        log(wrk)
    }else{
        wrk
    }
}

psi.Dn <- function(x, cdf = pnorm, type = c("D!=","D^+","D^-"),
                                   ties.jitter = NULL, ...){
    # if(missing(type))
    #     print("type is missing.")
    # else
    #     print(type)

    type <- match.arg(type)

    if (is.character(cdf))
        cdf <- get(cdf, mode = "function")
    if (mode(cdf) != "function")
        stop("'cdf' must be a function or a string naming a valid function")

    # podrobno razglezhdane na vaprosa e iztochnik na  chubavi zadachi. ???
    #
    # when jitter is done repeated calls of psi.Dn on the same vector will give slightly
    # different results.  ???
    #
    if(!is.null(ties.jitter)){
        if(isTRUE(ties.jitter))
            ties.jitter <- 0.5          # appropriate for data rounded to integer.
                                        # see jitter for other possible default
                                        # do  some  calcs to  determine ...
        if(!is.numeric(ties.jitter))
            stop("argument ties.jitter must be numeric or TRUE.")

        x <- psi.jitter(x,amount=0.5)  # modify x to get rid of ties.
    }


    n <- length(x)
    Fx <- cdf(sort(x),...)

    Dplus  <- max((1:n)/n - Fx)
    Dminus <- max(Fx  - (0:(n-1))/n)

    Dn <- max(Dplus,Dminus)

    switch(type, "D!=" = Dn,
                 "D^+" = Dplus,
                 "D^-" = Dminus)
}

# this assumes randgen and cdf accept the same parameters (supplied as ...).
psi.rks <- function(n,df,randgen = runif, cdf=punif, ...){
    # the commented out version givesan error, due  to using ... in psi.Dn
    # replicate(n, {x <- randgen(df,...); print(x); psi.Dn(x, cdf=cdf,...)} )
    # replicate(n, expression({x <- randgen(df,...); print(x); psi.Dn(x, cdf=cdf,...)}) )
    f <- function(){
        x <- randgen(df,...)
        psi.Dn(x, cdf=cdf)
    }
    replicate(n, f() )
}

psi.tied <- function(x){   # like duplicated but returns all tied elements
    x %in% unique(x[duplicated(x)])
}

psi.jitter <- function(x,amount){  # like jitter but jitters only the tied elements.
    x[psi.tied(x)] <- jitter(x[psi.tied(x)],amount=amount)
    x
}

psi.rlks.exp <- function(n,df){
    replicate(n, {x <- rexp(df); z <- x/mean(x); psi.Dn(z, pexp)})
}

psi.plks.exp <- function(q,df,Nsim=1000, lower.tail = TRUE){
    Dn  <- psi.rlks.exp(Nsim,df)
    FDn <- ecdf(Dn)                   # pomisli za zapazvane na rezultata za reuse. ???
    res <- FDn(q)                     # ??? obrabotka na  ties?
    res <- if(lower.tail) res else 1 - res
    names(res) <- paste(q)
    res
}

psi.qlks.exp <- function(p,df,Nsim=1000){
    Dn <- psi.rlks.exp(Nsim,df)
    res <- quantile(Dn, p)
    names(res) <- paste(p)
    res
}

psi.lks.exp.test <- function(x,Nsim=1000,...){
    n  <- length(x)
    z  <- x/mean(x)
    Dn <- psi.Dn(z, pexp,...)
    pval <- psi.plks.exp(Dn,n,Nsim=Nsim,lower.tail=FALSE,...)
    names(Dn) <- "Dn.Lillie.exp"

    # now make the return values similar to that of ks.test and other test functions.
    DNAME <- deparse(substitute(x))
    res <- list(statistic = Dn,
                p.value = pval,
                alternative = "two-sided",
                method = "One-sample Lilliefors test for exponential distribution",
                data.name = DNAME)
    class(res) <- "htest"
    res
}

# L for normal
# exp(-7.01256 * Kd^2 * (nd + 2.78019) + 2.99587 *
#         Kd * sqrt(nd + 2.78019) - 0.122119 + 0.974598/sqrt(nd) +
#         1.67997/nd)

## needs stepfun object  with finite domain.
##  val.right should really be NA.
psi.eqf <- function(x, val.right = NULL){
    x <- sort(x)
    if(is.null(val.right))
        val.right <- x[length(x)]
  # stepfun((1:(length(x)))/length(x), c(x, NA), right = TRUE)
    stepfun((1:(length(x)))/length(x), c(x, val.right), right = TRUE)
}

psi.plot.stepfun <-
function (x, xval, xlim, ylim = range(c(y, Fn.kn)), xlab = "x",
          ylab = "f(x)", main = NULL, add = FALSE, verticals = TRUE,
          do.points = TRUE, pch = par("pch"), col.points = par("col"),
          cex.points = par("cex"), col.hor = par("col"), col.vert = par("col"),
          lty = par("lty"), lwd = par("lwd"), rigid.xlim = FALSE, ...){
    if (!is.stepfun(x)) {
        if (is.numeric(x)) {
            sarg <- substitute(x)
            x <- ecdf(x)
            attr(x, "call") <- call("ecdf", sarg)
        }
        else stop("'plot.stepfun' called with wrong type of argument 'x'")
    }
    if (missing(main))
        main <- {
            cl <- attr(x, "call")
            deparse(if (!is.null(cl))
                cl
            else sys.call())
        }
    knF <- knots(x)
    xval <- if (missing(xval))
        knF
    else sort(xval)
    if (missing(xlim)) {
        rx <- range(xval)
        dr <- if (length(xval) > 1)
            max(0.08 * diff(rx), median(diff(xval)))
        else abs(xval)/16
        xlim <- rx + dr * c(-1, 1)
    }
    else dr <- diff(xlim)

    if(rigid.xlim){                                                             # Bosh ??? !!!
        knF <- knF[xlim[1] <= knF & knF <= xlim[2] ]
        ti <-  knF
        if(ti[1] != xlim[1]){
            ti <- c(xlim[1],ti)
        }
        if(ti[length(ti)] != xlim[2]){
            ti <- c(ti,xlim[2])
        }

    }else{
        knF <- knF[xlim[1] - dr <= knF & knF <= xlim[2] + dr]
        ti <- c(xlim[1] - dr, knF, xlim[2] + dr)
    }

    ti.l <- ti[-length(ti)]
    ti.r <- ti[-1]
    y <- x(0.5 * (ti.l + ti.r))
    n <- length(y)
    Fn.kn <- x(knF)
    if (add)
        segments(ti.l, y, ti.r, y, col = col.hor, lty = lty,
            lwd = lwd, ...)
    else {
        if (missing(ylim))
            ylim <- range(c(y, Fn.kn))
        plot(0, 0, type = "n", xlim = xlim, ylim = ylim, xlab = xlab,
            ylab = ylab, main = main, ...)
        segments(ti.l, y, ti.r, y, col = col.hor, lty = lty,
            lwd = lwd)
    }
    if (do.points)
        points(knF, Fn.kn, pch = pch, col = col.points, cex = cex.points)
    if (verticals){
        if(rigid.xlim){
            segments(ti.r[-n], y[-n], ti.r[-n], y[-1], col = col.vert,
                     lty = lty,
                     lwd = lwd)
        }else segments(knF, y[-n], knF, y[-1], col = col.vert, lty = lty,
                       lwd = lwd)
    }
    invisible(list(t = ti, y = y))
}