R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("accuracy-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('accuracy')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "PTB"
> 
> ### * PTB
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: PTBi
> ### Title: functions to add random noise to a vector
> ### Aliases: PTBi PTBms PTBmsb PTBmsbr PTBn PTBnc PTBns PTBnbr PTBnbrr
> ###   PTBnsbr PTBnsbrr PTBu PTBuc PTBus PTBubr PTBubrr PTBusbr PTBusbrr
> ### Keywords: misc optimize distribution
> 
> ### ** Examples
> 
>         x=1:1000
>         x.i = PTBi(x)
>         x.m = PTBms(x)
>         x.u = PTBu(x,size=.01)
>         x.us = PTBus(x,size=.01)
>         x.uc = PTBuc(x,size=.01)
> 
>         sum(x-x.i)              # should be 0, identity transform
[1] 0
>         sum(which(x!=x.i))      # also 0
[1] 0
> 
>         sum(x-x.m)              # should be quite small, very close to 0
[1] -1.464306e-11
>         sum(which(x!=x.m))      # should be near 1000
[1] 500500
> 
>         sum(x-x.u)              # should be small
[1] 0.0970383
>         sum(which(x!=x.m))      # should be near 1000
[1] 500500
> 
>         sum(x-x.us)             # should be small
[1] 36.09995
>         sum(which(x!=x.m))      # should be near 1000
[1] 500500
> 
>         sum(x-x.uc)             # should be near 0
[1] 1.196820e-13
>         sum(which(x!=x.m))      # should be near 1000
[1] 500500
>         
> 
> 
> 
> cleanEx(); ..nameEx <- "PTBm"
> 
> ### * PTBm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: PTBmu.gen
> ### Title: generator functions for multiple rounds of noise
> ### Aliases: PTBmu.gen PTBmn.gen
> ### Keywords: misc optimize distribution
> 
> ### ** Examples
> 
>         x=1:1000
>         f1=PTBmu.gen(); # should be roughly equivalent to PTBu() 
>         x.u = f1(x,size=1)
>         mean(x-x.u)  #should be small
[1] 0.0003083273
>         f2=PTBmu.gen(reps=100); # multiple disturbances tend to cancel eachother out
>         x.u2 = f2(x,size=1)
>         mean(x-x.u2)  #should be smaller 
[1] 0.0004959141
> 
> 
> 
> cleanEx(); ..nameEx <- "dehaan"
> 
> ### * dehaan
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: de Haan Global Optimality Test
> ### Title: Function to test that an optimum from MLE, NLS, or other
> ###   non-linear optization routine is a global optimum
> ### Aliases: dehaan
> ### Keywords: manip misc math htest
> 
> ### ** Examples
> 
> 
>   # The deHaan test is constructed as a maximum likelihood
>   # test, with negative values for the likelihood.  The BOD problem
>   # is a non-linear least squares minimization problem.  This test
>   # is implemented using the negative of the sum of squares for consistency
>   # with the deHaan framework of maximum likelihood.
> 
>   data("BOD",package="stats")
Warning in data("BOD", package = "stats") : datasets have been moved from package 'stats' to package 'datasets'
>   stval=expand.grid(A = seq(10, 100, 10), lrc = seq(.5, .8, .1))
>   llfun<-function(A,lrc)
+         -sum((BOD$demand - A*(1-exp(-exp(lrc)*BOD$Time)))^2)
>   lls=NULL 
>   for (i in 1:nrow(stval))  {
+      lls = rbind(lls, llfun(stval[i,1], stval[i,2]))
+   }
>   fm1 <- nls(demand ~ A*(1-exp(-exp(lrc)*Time)),  
+                 data = BOD, start = c(A = 20, lrc = log(.35)))
>   ss = -sum(resid(fm1)^2)
>   dehaan(lls, ss)
[1] TRUE
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "disttest"
> 
> ### * disttest
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Distribution Test Data, and log relative error function
> ### Title: Benchmark data to test the accuracy of statistical distribution
> ###   functions, function to compute log relative error
> ### Aliases: chisqtst ftst gammatst normtst ttst LRE
> ### Keywords: misc debugging
> 
> ### ** Examples
> 
> 
> # simple LRE examples
> LRE(1.001,1) # roughly 3 significant digits agreement
[1] 3
> LRE(1,1) # complete agreement
[1] Inf
> LRE(20,1) # complete disagreement
[1] -1.278754
> 
> #
> # how accurate are student's t-test functions?
> #
> 
> data(ttst)
> # compute t quantiles using benchmark data
> tqt = qt(ttst$p,ttst$df)
> 
> # compute log-relative-error (LRE) of  qt() results, compared to 
> # correct answers
> 
> lrq = LRE(tqt, ttst$invt);
Warning in log(x, base) : NaNs produced
> 
> # if there are entries with LRE's of < 5, there may be
> # significant inaccuracies in the qt() function
> 
> table(trunc(lrq))

-Inf    3    4    5    6    7    8    9  Inf 
   1    1   13  814 2536  322   18    4   34 
> 
> # now repeat process, for pt()
> 
> tpt = pt(ttst$invt,ttst$df)
> lrp= LRE(tpt, ttst$pinvt);
> table(trunc(lrp))

   5    6    7    8    9   10   11   12   13  Inf 
 636 3820 1071  556  492  459  358   23    1   36 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "frexp"
> 
> ### * frexp
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: frexp
> ### Title: Function to convert vector of floating-point numbers to
> ###   fractional and integral components
> ### Aliases: frexp
> ### Keywords: manip misc math
> 
> ### ** Examples
> 
> x = runif(10)
> y = frexp(x)
> y
       Mantissa Exponent
 [1,] 0.5310173       -1
 [2,] 0.7442478       -1
 [3,] 0.5728534        0
 [4,] 0.9082078        0
 [5,] 0.8067277       -2
 [6,] 0.8983897        0
 [7,] 0.9446753        0
 [8,] 0.6607978        0
 [9,] 0.6291140        0
[10,] 0.9885803       -4
> # this is now true by construction
> x==y[,1] *2^y[,2]
 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> 
> 
> 
> cleanEx(); ..nameEx <- "perturb"
> 
> ### * perturb
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: perturb
> ### Title: Data Perturbations based Sensitivity Analysis
> ### Aliases: perturb
> ### Keywords: misc optimize
> 
> ### ** Examples
> 
> 
> # Examine the sensitivity of the GLM from Venables & Ripley (2002, p.189)
> # as described in the glm module.
> # 
> # Perturb the two independent variables using +/- 0.025
> #       (relative to the size of each observations)
> # uniformly distributed noise. Dependent variable is not being modified.
> # 
> # Summary should show that estimated coefficients are not substantively affected by noise.
> 
> data(anorexia,package="MASS")
> panorexia = perturb(anorexia, glm, Postwt ~ Prewt + Treat + offset(Prewt),
+      family=gaussian, 
+     ptb.R=100, ptb.ran.gen=c(PTBi,PTBus,PTBus), ptb.s=c(1,.005,.005) )
> summary(panorexia)
statistic: 
function (formula, family = gaussian, data, weights, subset, 
    na.action, start = NULL, etastart, mustart, offset, control = glm.control(...), 
    model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL, 
    ...) 
{
    call <- match.call()
    if (is.character(family)) 
        family <- get(family, mode = "function", envir = parent.frame())
    if (is.function(family)) 
        family <- family()
    if (is.null(family$family)) {
        print(family)
        stop("'family' not recognized")
    }
    if (missing(data)) 
        data <- environment(formula)
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data", "subset", "weights", "na.action", 
        "etastart", "mustart", "offset"), names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
    switch(method, model.frame = return(mf), glm.fit = 1, glm.fit.null = 1, 
        stop("invalid 'method': ", method))
    mt <- attr(mf, "terms")
    Y <- model.response(mf, "numeric")
    if (length(dim(Y)) == 1) {
        nm <- rownames(Y)
        dim(Y) <- NULL
        if (!is.null(nm)) 
            names(Y) <- nm
    }
    X <- if (!is.empty.model(mt)) 
        model.matrix(mt, mf, contrasts)
    else matrix(, NROW(Y), 0)
    weights <- model.weights(mf)
    offset <- model.offset(mf)
    if (!is.null(weights) && any(weights < 0)) 
        stop("negative weights not allowed")
    if (!is.null(offset) && length(offset) != NROW(Y)) 
        stop(gettextf("number of offsets is %d should equal %d (number of observations)", 
            length(offset), NROW(Y)), domain = NA)
    mustart <- model.extract(mf, "mustart")
    etastart <- model.extract(mf, "etastart")
    fit <- glm.fit(x = X, y = Y, weights = weights, start = start, 
        etastart = etastart, mustart = mustart, offset = offset, 
        family = family, control = control, intercept = attr(mt, 
            "intercept") > 0)
    if (any(offset) && attr(mt, "intercept") > 0) {
        fit$null.deviance <- glm.fit(x = X[, "(Intercept)", drop = FALSE], 
            y = Y, weights = weights, offset = offset, family = family, 
            control = control, intercept = TRUE)$deviance
    }
    if (model) 
        fit$model <- mf
    fit$na.action <- attr(mf, "na.action")
    if (x) 
        fit$x <- X
    if (!y) 
        fit$y <- NULL
    fit <- c(fit, list(call = call, formula = formula, terms = mt, 
        data = data, offset = offset, control = control, method = method, 
        contrasts = attr(X, "contrasts"), xlevels = .getXlevels(mt, 
            mf)))
    class(fit) <- c("glm", "lm")
    fit
}
<environment: namespace:stats>
s: 
[1] 1.000 0.005 0.005
Replications: 
 100 

ran.gen: 
[[1]]
function (x, size = 1) 
{
    return(x)
}
<environment: namespace:accuracy>

[[2]]
function (x, size = 1) 
{
    return(PTBunif(x, size = size, scaled = TRUE))
}
<environment: namespace:accuracy>

[[3]]
function (x, size = 1) 
{
    return(PTBunif(x, size = size, scaled = TRUE))
}
<environment: namespace:accuracy>

formula: 
 

Postwt ~ Prewt + Treat + offset(Prewt)
attr(,"variables")
list(Postwt, Prewt, Treat, offset(Prewt))
attr(,"offset")
[1] 4
attr(,"factors")
              Prewt Treat
Postwt            0     0
Prewt             1     0
Treat             0     1
offset(Prewt)     0     0
attr(,"term.labels")
[1] "Prewt" "Treat"
attr(,"order")
[1] 1 1
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: R_GlobalEnv>
attr(,"predvars")
list(Postwt, Prewt, Treat, offset(Prewt))
attr(,"dataClasses")
       Postwt         Prewt         Treat offset(Prewt) 
    "numeric"     "numeric"      "factor"     "numeric" 
betas:

  X.Intercept.       Prewt           TreatCont         TreatFT     
 Min.   :48.86   Min.   :-0.5773   Min.   :-4.208   Min.   :4.459  
 1st Qu.:49.62   1st Qu.:-0.5689   1st Qu.:-4.131   1st Qu.:4.526  
 Median :49.86   Median :-0.5667   Median :-4.099   Median :4.556  
 Mean   :49.82   Mean   :-0.5661   Mean   :-4.098   Mean   :4.557  
 3rd Qu.:50.05   3rd Qu.:-0.5634   3rd Qu.:-4.071   3rd Qu.:4.587  
 Max.   :50.70   Max.   :-0.5546   Max.   :-3.977   Max.   :4.654  
      (Intercept)      Prewt TreatCont  TreatFT
2.5%     48.97372 -0.5749067 -4.187017 4.465082
97.5%    50.53673 -0.5555569 -4.010346 4.641438
stderrs:

  X.Intercept.       Prewt          TreatCont        TreatFT     
 Min.   :13.18   Min.   :0.1586   Min.   :1.884   Min.   :2.123  
 1st Qu.:13.34   1st Qu.:0.1606   1st Qu.:1.890   1st Qu.:2.130  
 Median :13.39   Median :0.1612   Median :1.894   Median :2.134  
 Mean   :13.39   Mean   :0.1611   Mean   :1.894   Mean   :2.134  
 3rd Qu.:13.42   3rd Qu.:0.1616   3rd Qu.:1.896   3rd Qu.:2.136  
 Max.   :13.53   Max.   :0.1628   Max.   :1.906   Max.   :2.146  
      (Intercept)     Prewt TreatCont  TreatFT
2.5%     13.27673 0.1597773  1.886697 2.125414
97.5%    13.48918 0.1624003  1.902758 2.142993
> 
> # Use classic longley dataset. The model is numerically unstable, 
> # and much more sensitive to noise.  Smaller amounts of noise tremendously
> # alter some of the estimated coefficients:
> # 
> # In this example we are not perturbing the dependent variable (employed) or 
> # the year variable. So we assign then PTBi or NULL in ptb.ran.gen )
> 
> if(version$major<2) {
+         data(longley, package="base")
+ } else {
+         data(longley, package="datasets")
+ }
> 
> plongley = perturb(longley,lm,Employed~., ptb.R=100, 
+     ptb.ran.gen=c(PTBi, replicate(5,PTBus),PTBi), ptb.s=c(1,replicate(5,.001),1))
> sp=summary(plongley)
> 
> # summarizes range
> print(sp)
statistic: 
function (formula, data, subset, weights, na.action, method = "qr", 
    model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
    contrasts = NULL, offset, ...) 
{
    ret.x <- x
    ret.y <- y
    cl <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data", "subset", "weights", "na.action", 
        "offset"), names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
    if (method == "model.frame") 
        return(mf)
    else if (method != "qr") 
        warning(gettextf("method = '%s' is not supported. Using 'qr'", 
            method), domain = NA)
    mt <- attr(mf, "terms")
    y <- model.response(mf, "numeric")
    w <- model.weights(mf)
    offset <- model.offset(mf)
    if (!is.null(offset) && length(offset) != NROW(y)) 
        stop(gettextf("number of offsets is %d, should equal %d (number of observations)", 
            length(offset), NROW(y)), domain = NA)
    if (is.empty.model(mt)) {
        x <- NULL
        z <- list(coefficients = if (is.matrix(y)) matrix(, 0, 
            3) else numeric(0), residuals = y, fitted.values = 0 * 
            y, weights = w, rank = 0, df.residual = if (is.matrix(y)) nrow(y) else length(y))
        if (!is.null(offset)) 
            z$fitted.values <- offset
    }
    else {
        x <- model.matrix(mt, mf, contrasts)
        z <- if (is.null(w)) 
            lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)
        else lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, 
            ...)
    }
    class(z) <- c(if (is.matrix(y)) "mlm", "lm")
    z$na.action <- attr(mf, "na.action")
    z$offset <- offset
    z$contrasts <- attr(x, "contrasts")
    z$xlevels <- .getXlevels(mt, mf)
    z$call <- cl
    z$terms <- mt
    if (model) 
        z$model <- mf
    if (ret.x) 
        z$x <- x
    if (ret.y) 
        z$y <- y
    if (!qr) 
        z$qr <- NULL
    z
}
<environment: namespace:stats>
s: 
[1] 1.000 0.001 0.001 0.001 0.001 0.001 1.000
Replications: 
 100 

ran.gen: 
[[1]]
function (x, size = 1) 
{
    return(x)
}
<environment: namespace:accuracy>

[[2]]
function (x, size = 1) 
{
    return(PTBunif(x, size = size, scaled = TRUE))
}
<environment: namespace:accuracy>

[[3]]
function (x, size = 1) 
{
    return(PTBunif(x, size = size, scaled = TRUE))
}
<environment: namespace:accuracy>

[[4]]
function (x, size = 1) 
{
    return(PTBunif(x, size = size, scaled = TRUE))
}
<environment: namespace:accuracy>

[[5]]
function (x, size = 1) 
{
    return(PTBunif(x, size = size, scaled = TRUE))
}
<environment: namespace:accuracy>

[[6]]
function (x, size = 1) 
{
    return(PTBunif(x, size = size, scaled = TRUE))
}
<environment: namespace:accuracy>

[[7]]
function (x, size = 1) 
{
    return(x)
}
<environment: namespace:accuracy>

formula: 
 

statistic(formula = ..1, data = ndata)
betas:

  X.Intercept.       GNP.deflator           GNP             Unemployed       
 Min.   :-1539.01   Min.   :-0.17319   Min.   :-0.05813   Min.   :-0.022694  
 1st Qu.: -622.24   1st Qu.:-0.06098   1st Qu.: 0.05205   1st Qu.:-0.006999  
 Median : -317.20   Median :-0.04789   Median : 0.06258   Median :-0.005445  
 Mean   : -302.42   Mean   :-0.04345   Mean   : 0.06054   Mean   :-0.005846  
 3rd Qu.:   12.12   3rd Qu.:-0.02532   3rd Qu.: 0.07003   3rd Qu.:-0.004279  
 Max.   :  993.84   Max.   : 0.30254   Max.   : 0.13001   Max.   : 0.002829  
  Armed.Forces         Population           Year         
 Min.   :-0.011682   Min.   :-0.8457   Min.   :-0.46865  
 1st Qu.:-0.006564   1st Qu.:-0.4225   1st Qu.: 0.04074  
 Median :-0.005919   Median :-0.3836   Median : 0.20774  
 Mean   :-0.006139   Mean   :-0.3669   Mean   : 0.20222  
 3rd Qu.:-0.005506   3rd Qu.:-0.2954   3rd Qu.: 0.36623  
 Max.   :-0.002024   Max.   : 0.5789   Max.   : 0.83550  
      (Intercept) GNP.deflator        GNP    Unemployed Armed.Forces Population
2.5%   -1375.9575  -0.14211591 0.02053138 -0.0113143147 -0.008764543 -0.6856307
97.5%    847.7891   0.09246326 0.10457074 -0.0002767496 -0.004649820 -0.1039797
            Year
2.5%  -0.3780163
97.5%  0.7463708
stderrs:

  X.Intercept.    GNP.deflator          GNP            Unemployed      
 Min.   :243.2   Min.   :0.08316   Min.   :0.02054   Min.   :0.002787  
 1st Qu.:404.2   1st Qu.:0.12417   1st Qu.:0.03101   1st Qu.:0.004394  
 Median :493.1   Median :0.13741   Median :0.03403   Median :0.004778  
 Mean   :492.8   Mean   :0.13415   Mean   :0.03552   Mean   :0.005006  
 3rd Qu.:567.4   3rd Qu.:0.14425   3rd Qu.:0.03828   3rd Qu.:0.005498  
 Max.   :864.7   Max.   :0.18737   Max.   :0.05802   Max.   :0.008071  
  Armed.Forces        Population          Year       
 Min.   :0.001804   Min.   :0.2037   Min.   :0.1250  
 1st Qu.:0.002668   1st Qu.:0.3105   1st Qu.:0.2083  
 Median :0.002949   Median :0.3430   Median :0.2480  
 Mean   :0.002919   Mean   :0.3395   Mean   :0.2509  
 3rd Qu.:0.003137   3rd Qu.:0.3587   3rd Qu.:0.2887  
 Max.   :0.004222   Max.   :0.4944   Max.   :0.4404  
      (Intercept) GNP.deflator        GNP  Unemployed Armed.Forces Population
2.5%     281.6912    0.0877617 0.02244478 0.003234530  0.001992270  0.2267091
97.5%    683.3782    0.1669729 0.05166461 0.007281459  0.003764202  0.4597971
           Year
2.5%  0.1430853
97.5% 0.3514665
> 
> # look in summary object to extract more ...
> names(attributes(sp))
[1] "coef.names.m"   "coef.betas.m"   "coef.stderrs.m" "coef.formula"  
[5] "ran.gen"        "R"              "s"              "statistic"     
[9] "class"         
> 
> # print metrix of coefficients from all runs
> coef= attr(sp,"coef.betas.m")
> summary(coef)
  X.Intercept.       GNP.deflator           GNP             Unemployed       
 Min.   :-1539.01   Min.   :-0.17319   Min.   :-0.05813   Min.   :-0.022694  
 1st Qu.: -622.24   1st Qu.:-0.06098   1st Qu.: 0.05205   1st Qu.:-0.006999  
 Median : -317.20   Median :-0.04789   Median : 0.06258   Median :-0.005445  
 Mean   : -302.42   Mean   :-0.04345   Mean   : 0.06054   Mean   :-0.005846  
 3rd Qu.:   12.12   3rd Qu.:-0.02532   3rd Qu.: 0.07003   3rd Qu.:-0.004279  
 Max.   :  993.84   Max.   : 0.30254   Max.   : 0.13001   Max.   : 0.002829  
  Armed.Forces         Population           Year         
 Min.   :-0.011682   Min.   :-0.8457   Min.   :-0.46865  
 1st Qu.:-0.006564   1st Qu.:-0.4225   1st Qu.: 0.04074  
 Median :-0.005919   Median :-0.3836   Median : 0.20774  
 Mean   :-0.006139   Mean   :-0.3669   Mean   : 0.20222  
 3rd Qu.:-0.005506   3rd Qu.:-0.2954   3rd Qu.: 0.36623  
 Max.   :-0.002024   Max.   : 0.5789   Max.   : 0.83550  
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "sechol"
> 
> ### * sechol
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sechol
> ### Title: Schnabel-Eskow Choleksy Decomposition
> ### Aliases: sechol
> ### Keywords: array optimize
> 
> ### ** Examples
> 
> # compare with chol() on a non-singular matrix
> S <- matrix(c(2,0,2.4,0,2,0,2.4,0,3),ncol=3)
> chol(S)
         [,1]     [,2]      [,3]
[1,] 1.414214 0.000000 1.6970563
[2,] 0.000000 1.414214 0.0000000
[3,] 0.000000 0.000000 0.3464102
> T <- sechol(S)
> t(T) 
         [,1]     [,2]      [,3]
[1,] 1.414214 0.000000 0.0000000
[2,] 0.000000 1.414214 0.0000000
[3,] 1.697056 0.000000 0.3464102
attr(,"delta")
[1] 0
> 
> # an example with a singular matrix
> S <- matrix(c(2,0,2.5,0,2,0,2.5,0,3),ncol=3)
> sechol(S) 
         [,1]    [,2]        [,3]
[1,] 1.414214 0.00000 1.767766953
[2,] 0.000000 1.41422 0.000000000
[3,] 0.000000 0.00000 0.004262202
attr(,"delta")
[1] 1.816636e-05
> t(T) 
         [,1]     [,2]      [,3]
[1,] 1.414214 0.000000 0.0000000
[2,] 0.000000 1.414214 0.0000000
[3,] 1.697056 0.000000 0.3464102
attr(,"delta")
[1] 0
>     
> # another example with a singular matrix
> S <- matrix(c(2,0,10,0,2,0,10,0,3),ncol=3)
> T <- sechol(S)
> t(T) 
         [,1]     [,2]       [,3]
[1,] 1.414214 0.000000 0.00000000
[2,] 0.000000 1.414235 0.00000000
[3,] 7.071068 0.000000 0.00778168
attr(,"delta")
[1] 6.055454e-05
> 
> 
> 
> cleanEx(); ..nameEx <- "starr"
> 
> ### * starr
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Starr Global Optimality Test
> ### Title: Function to test that an optimum from MLE, NLS, or other
> ###   non-linear optization routine is a global optimum
> ### Aliases: starr
> ### Keywords: manip misc math htest
> 
> ### ** Examples
> 
> 
>         x=rbind(c(1,1,1), c(1,2,1), c(1,1.1,1), c(1,2,1), c(3,4,5));
>         starr(rbind(1,1,2,2));
[1] 0.3333333
> 
> 
> cleanEx(); ..nameEx <- "trueRandom"
> 
> ### * trueRandom
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: trueRandom
> ### Title: Function to return TRUE (not pseudo) random numbers, based on
> ###   system and networked entropy collection.
> ### Aliases: trueRandom resetSeed runifT initPool
> ### Keywords: misc datagen distribution
> 
> ### ** Examples
> 
> # note, if you are using Windows, you must be on-line
> # to get to entropy generator, will fall-back to pseudo-random
> # numbers when off-line
> 
> # only needed  to pass R CMD check on online windows
> w = options("warn");
> options(warn=0);
> 
> # reset the seed for runif() based on a true random value
> 
> resetSeed()
> y=runif(100)
> ty= table(trunc(5*y))
> print(ty);

 0  1  2  3  4 
25 19 20 13 23 
> chisq.test(ty)

	Chi-squared test for given probabilities

data:  ty 
X-squared = 4.2, df = 4, p-value = 0.3796

> 
> # generate true random values directly (may block for long periods if 
> # if entropy pool is empty)
> 
> y=runifT(100)
> ty= table(trunc(5*y))
> print(ty);

 0  1  2  3  4 
21 21 20 16 22 
> chisq.test(ty)

	Chi-squared test for given probabilities

data:  ty 
X-squared = 1.1, df = 4, p-value = 0.8943

> options(warn=as.integer(w));
> 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  11.53 0.17 12.1 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')