Full_Name: Madeleine Yeh
Version: 1.9.1
OS: AIX 5.2
Submission from: (NULL) (151.121.225.1)
After compiling R-1.9.1 on AIX 5.2 using the IBM cc compiler, I ran the
checks. One of them failed. Here is the output from running the check solo.
root@svweb:/fsapps/test/build/R/1.9.1/R-1.9.1/tests/Examples:># ../../bin/R --vanilla < stats-Ex.R
R : Copyright 2004, The R Foundation for Statistical Computing
Version 1.9.1 (2004-06-21), ISBN 3-900051-00-3
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 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, at least for base graphics
> .newplot.hook <- function()
+ {
+ pp <-
par(c("mfg","mfcol","oma","mar"))
+ if(all(pp$mfg[1:2] == c(1, pp$mfcol[2]))) {
+ outer <- (oma4 <- pp$oma[4]) > 0; mar4 <- pp$mar[4]
+ mtext(paste("help(", ..nameEx, ")"), side = 4,
+ line = if(outer)max(1, oma4 - 1) else min(1, mar4 - 1),
+ outer = outer, adj = 1, cex = .8, col = "orchid", las=3)
+ }
+ }> setHook("plot.new", .newplot.hook)
> setHook("persp", .newplot.hook)
> rm(.newplot.hook)
> ## add some hooks to label plot pages for grid graphics
> .gridplot.hook <- function()
+ {
+ pushViewport(viewport(width=unit(1, "npc") - unit(1,
"lines"),
+ x=0, just="left"))
+ grid.text(paste("help(", ..nameEx, ")"),
+ x=unit(1, "npc") + unit(0.5, "lines"),
+ y=unit(0.8, "npc"), rot=90,
+ gp=gpar(col="orchid"))
+ }> setHook("grid.newpage", .gridplot.hook)
> rm(.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)
+ assign("T", delay(stop("T used instead of
TRUE")),
+ pos = .CheckExEnv)
+ assign("F", delay(stop("F used instead of
FALSE")),
+ pos = .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)
> graphics::postscript("stats-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env
.CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered =
"contr.poly"))
> library('stats')
>
> assign(".oldSearch", search(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "AIC"
>
> ### * AIC
>
> flush(stderr()); flush(stdout())
>
> ### Name: AIC
> ### Title: Akaike's An Information Criterion
> ### Aliases: AIC
> ### Keywords: models
>
> ### ** Examples
>
> data(swiss)
> lm1 <- lm(Fertility ~ . , data = swiss)
> AIC(lm1)
[1] 326.0716> stopifnot(all.equal(AIC(lm1),
+ AIC(logLik(lm1))))> ## a version of BIC or Schwarz' BC :
> AIC(lm1, k = log(nrow(swiss)))
[1] 339.0226>
>
>
> cleanEx(); ..nameEx <- "ARMAacf"
>
> ### * ARMAacf
>
> flush(stderr()); flush(stdout())
>
> ### Name: ARMAacf
> ### Title: Compute Theoretical ACF for an ARMA Process
> ### Aliases: ARMAacf
> ### Keywords: ts
>
> ### ** Examples
>
> ARMAacf(c(1.0, -0.25), 1.0, lag.max = 10)
0 1 2 3 4 5
1.000000000 0.875000000 0.625000000 0.406250000 0.250000000 0.148437500
6 7 8 9 10
0.085937500 0.048828125 0.027343750 0.015136719 0.008300781
> ## Example from Brockwell & Davis (1991, pp.92-4)
> ## answer 2^(-n) * (32/3 + 8 * n) /(32/3)
> n <- 1:10; 2^(-n) * (32/3 + 8 * n) /(32/3)
[1] 0.875000000 0.625000000 0.406250000 0.250000000 0.148437500 0.085937500
[7] 0.048828125 0.027343750 0.015136719 0.008300781> ARMAacf(c(1.0, -0.25), 1.0, lag.max = 10, pacf = TRUE)
[1] 0.8750000 -0.6000000 0.3750000 -0.2727273 0.2142857 -0.1764706
[7] 0.1500000 -0.1304348 0.1153846 -0.1034483> ARMAacf(c(1.0, -0.25), lag.max = 10, pacf = TRUE)
[1] 8.000000e-01 -2.500000e-01 -6.579099e-17 -3.045879e-19 5.302086e-19
[6] 2.059071e-17 -6.361500e-33 -1.027984e-17 5.139921e-18
1.606263e-33>
>
>
> cleanEx(); ..nameEx <- "ARMAtoMA"
>
> ### * ARMAtoMA
>
> flush(stderr()); flush(stdout())
>
> ### Name: ARMAtoMA
> ### Title: Convert ARMA Process to Infinite MA Process
> ### Aliases: ARMAtoMA
> ### Keywords: ts
>
> ### ** Examples
>
> ARMAtoMA(c(1.0, -0.25), 1.0, 10)
[1] 2.00000000 1.75000000 1.25000000 0.81250000 0.50000000 0.29687500
[7] 0.17187500 0.09765625 0.05468750 0.03027344> ## Example from Brockwell & Davis (1991, p.92)
> ## answer (1 + 3*n)*2^(-n)
> n <- 1:10; (1 + 3*n)*2^(-n)
[1] 2.00000000 1.75000000 1.25000000 0.81250000 0.50000000 0.29687500
[7] 0.17187500 0.09765625 0.05468750 0.03027344>
>
>
> cleanEx(); ..nameEx <- "AirPassengers"
>
> ### * AirPassengers
>
> flush(stderr()); flush(stdout())
>
> ### Name: AirPassengers
> ### Title: Monthly Airline Passenger Numbers 1949-1960
> ### Aliases: AirPassengers
> ### Keywords: datasets
>
> ### ** Examples
>
> ## Not run: ## These are quite slow and so not run by
example(AirPassengers)
> ##D
> ##D data(AirPassengers)
> ##D ## The classic 'airline model', by full ML
> ##D (fit <- arima(log10(AirPassengers), c(0, 1, 1),
> ##D seasonal = list(order=c(0, 1 ,1), period=12)))
> ##D update(fit, method = "CSS")
> ##D update(fit, x=window(log10(AirPassengers), start = 1954))
> ##D pred <- predict(fit, n.ahead = 24)
> ##D tl <- pred$pred - 1.96 * pred$se
> ##D tu <- pred$pred + 1.96 * pred$se
> ##D ts.plot(AirPassengers, 10^tl, 10^tu, log = "y", lty =
c(1,2,2))
> ##D
> ##D ## full ML fit is the same if the series is reversed, CSS fit is not
> ##D ap0 <- rev(log10(AirPassengers))
> ##D attributes(ap0) <- attributes(AirPassengers)
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order=c(0, 1 ,1), period=12))
> ##D arima(ap0, c(0, 1, 1), seasonal = list(order=c(0, 1 ,1), period=12),
> ##D method = "CSS")
> ##D
> ##D ## Structural Time Series
> ##D ap <- log10(AirPassengers) - 2
> ##D (fit <- StructTS(ap, type= "BSM"))
> ##D par(mfrow=c(1,2))
> ##D plot(cbind(ap, fitted(fit)), plot.type = "single")
> ##D plot(cbind(ap, tsSmooth(fit)), plot.type = "single")
> ## End(Not run)
>
>
> cleanEx(); ..nameEx <- "BOD"
>
> ### * BOD
>
> flush(stderr()); flush(stdout())
>
> ### Name: BOD
> ### Title: Biochemical Oxygen Demand
> ### Aliases: BOD
> ### Keywords: datasets
>
> ### ** Examples
>
> data(BOD)
> # simplest form of fitting a first-order model to these data
> fm1 <- nls(demand ~ A*(1-exp(-exp(lrc)*Time)), data = BOD,
+ start = c(A = 20, lrc = log(.35)))> coef(fm1)
A lrc
19.1425768 -0.6328215 > print(fm1)
Nonlinear regression model
model: demand ~ A * (1 - exp(-exp(lrc) * Time))
data: BOD
A lrc
19.1425768 -0.6328215
residual sum-of-squares: 25.99027 > # using the plinear algorithm
> fm2 <- nls(demand ~ (1-exp(-exp(lrc)*Time)), data = BOD,
+ start = c(lrc = log(.35)), algorithm = "plinear", trace = TRUE)
32.94622 : -1.049822 22.126001
25.99248 : -0.6257161 19.1031883
25.99027 : -0.6327039 19.1419223
25.99027 : -0.6328192 19.1425644 > # using a self-starting model
> fm3 <- nls(demand ~ SSasympOrig(Time, A, lrc), data = BOD)
> summary( fm3 )
Formula: demand ~ SSasympOrig(Time, A, lrc)
Parameters:
Estimate Std. Error t value Pr(>|t|)
A 19.1426 2.4959 7.670 0.00155 **
lrc -0.6328 0.3824 -1.655 0.17328
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
Residual standard error: 2.549 on 4 degrees of freedom
Correlation of Parameter Estimates:
A
lrc -0.8528
>
>
>
> cleanEx(); ..nameEx <- "Beta"
>
> ### * Beta
>
> flush(stderr()); flush(stdout())
>
> ### Name: Beta
> ### Title: The Beta Distribution
> ### Aliases: Beta dbeta pbeta qbeta rbeta
> ### Keywords: distribution
>
> ### ** Examples
>
> x <- seq(0, 1, length=21)
> dbeta(x, 1, 1)
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1> pbeta(x, 1, 1)
[1] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70
[16] 0.75 0.80 0.85 0.90 0.95 1.00>
>
>
> cleanEx(); ..nameEx <- "Binomial"
>
> ### * Binomial
>
> flush(stderr()); flush(stdout())
>
> ### Name: Binomial
> ### Title: The Binomial Distribution
> ### Aliases: Binomial dbinom pbinom qbinom rbinom
> ### Keywords: distribution
>
> ### ** Examples
>
> # Compute P(45 < X < 55) for X Binomial(100,0.5)
> sum(dbinom(46:54, 100, 0.5))
[1] 0.6317984>
> ## Using "log = TRUE" for an extended range :
> n <- 2000
> k <- seq(0, n, by = 20)
> plot (k, dbinom(k, n, pi/10, log=TRUE), type='l', ylab="log
density",
+ main = "dbinom(*, log=TRUE) is better than
log(dbinom(*))")> lines(k, log(dbinom(k, n, pi/10)), col='red', lwd=2)
> ## extreme points are omitted since dbinom gives 0.
> mtext("dbinom(k, log=TRUE)", adj=0)
> mtext("extended range", adj=0, line = -1, font=4)
> mtext("log(dbinom(k))", col="red", adj=1)
>
>
>
> cleanEx(); ..nameEx <- "Cauchy"
>
> ### * Cauchy
>
> flush(stderr()); flush(stdout())
>
> ### Name: Cauchy
> ### Title: The Cauchy Distribution
> ### Aliases: Cauchy dcauchy pcauchy qcauchy rcauchy
> ### Keywords: distribution
>
> ### ** Examples
>
> dcauchy(-1:4)
[1] 0.15915494 0.31830989 0.15915494 0.06366198 0.03183099
0.01872411>
>
>
> cleanEx(); ..nameEx <- "ChickWeight"
>
> ### * ChickWeight
>
> flush(stderr()); flush(stdout())
>
> ### Name: ChickWeight
> ### Title: Weight versus age of chicks on different diets
> ### Aliases: ChickWeight
> ### Keywords: datasets
>
> ### ** Examples
>
> data(ChickWeight)
> coplot(weight ~ Time | Chick, data = ChickWeight,
+ type = "b", show = FALSE)> ## fit a representative chick
> fm1 <- nls(weight ~ SSlogis( Time, Asym, xmid, scal ),
+ data = ChickWeight, subset = Chick == 1)> summary( fm1 )
Formula: weight ~ SSlogis(Time, Asym, xmid, scal)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Asym 937.0213 465.8578 2.011 0.07516 .
xmid 35.2228 8.3119 4.238 0.00218 **
scal 11.4052 0.9052 12.599 5.08e-07 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
Residual standard error: 2.919 on 9 degrees of freedom
Correlation of Parameter Estimates:
Asym xmid
xmid 0.9991
scal 0.9745 0.9829
>
>
>
> cleanEx(); ..nameEx <- "Chisquare"
>
> ### * Chisquare
>
> flush(stderr()); flush(stdout())
>
> ### Name: Chisquare
> ### Title: The (non-central) Chi-Squared Distribution
> ### Aliases: Chisquare dchisq pchisq qchisq rchisq
> ### Keywords: distribution
>
> ### ** Examples
>
> dchisq(1, df=1:3)
[1] 0.2419707 0.3032653 0.2419707> pchisq(1, df= 3)
[1] 0.1987480> pchisq(1, df= 3, ncp = 0:4)# includes the above
[1] 0.19874804 0.13229855 0.08787311 0.05824691
0.03853592>
> x <- 1:10
> ## Chi-squared(df = 2) is a special exponential distribution
> all.equal(dchisq(x, df=2), dexp(x, 1/2))
[1] TRUE> all.equal(pchisq(x, df=2), pexp(x, 1/2))
[1] TRUE>
> ## non-central RNG -- df=0 is ok for ncp > 0: Z0 has point mass at 0!
> Z0 <- rchisq(100, df = 0, ncp = 2.)
> graphics::stem(Z0)
The decimal point is at the |
0 | 0000000000000000000000000000000000000013356778899
1 | 0001333456678888899
2 | 0011444467
3 | 00233345888
4 | 111246
5 |
6 |
7 | 178
8 | 23
>
> ## Not run:
> ##D ## visual testing
> ##D ## do P-P plots for 1000 points at various degrees of freedom
> ##D L <- 1.2; n <- 1000; pp <- ppoints(n)
> ##D op <- par(mfrow = c(3,3), mar= c(3,3,1,1)+.1, mgp= c(1.5,.6,0),
> ##D oma = c(0,0,3,0))
> ##D for(df in 2^(4*rnorm(9))) {
> ##D plot(pp, sort(pchisq(rr <- rchisq(n,df=df, ncp=L), df=df, ncp=L)),
> ##D ylab="pchisq(rchisq(.),.)", pch=".")
> ##D mtext(paste("df = ",formatC(df, digits = 4)), line= -2,
adj=0.05)
> ##D abline(0,1,col=2)
> ##D }
> ##D mtext(expression("P-P plots : Noncentral "*
> ##D chi^2 *"(n=1000, df=X, ncp= 1.2)"),
> ##D cex = 1.5, font = 2, outer=TRUE)
> ##D par(op)
> ## End(Not run)
>
>
>
> cleanEx(); ..nameEx <- "DNase"
>
> ### * DNase
>
> flush(stderr()); flush(stdout())
>
> ### Name: DNase
> ### Title: Elisa assay of DNase
> ### Aliases: DNase
> ### Keywords: datasets
>
> ### ** Examples
>
> data(DNase)
> coplot(density ~ conc | Run, data = DNase,
+ show = FALSE, type = "b")> coplot(density ~ log(conc) | Run, data = DNase,
+ show = FALSE, type = "b")> ## fit a representative run
> fm1 <- nls(density ~ SSlogis( log(conc), Asym, xmid, scal ),
+ data = DNase, subset = Run == 1)> ## compare with a four-parameter logistic
> fm2 <- nls(density ~ SSfpl( log(conc), A, B, xmid, scal ),
+ data = DNase, subset = Run == 1)> summary(fm2)
Formula: density ~ SSfpl(log(conc), A, B, xmid, scal)
Parameters:
Estimate Std. Error t value Pr(>|t|)
A -0.007897 0.017200 -0.459 0.654
B 2.377239 0.109516 21.707 5.35e-11 ***
xmid 1.507403 0.102080 14.767 4.65e-09 ***
scal 1.062579 0.056996 18.643 3.16e-10 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
Residual standard error: 0.01981 on 12 degrees of freedom
Correlation of Parameter Estimates:
A B xmid
B -0.6375
xmid -0.5176 0.9811
scal -0.8030 0.9266 0.8796
> anova(fm1, fm2)
Analysis of Variance Table
Model 1: density ~ SSlogis(log(conc), Asym, xmid, scal)
Model 2: density ~ SSfpl(log(conc), A, B, xmid, scal)
Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1 13 0.0047896
2 12 0.0047073 1 0.0000823 0.2098 0.6551>
>
>
> cleanEx(); ..nameEx <- "Exponential"
>
> ### * Exponential
>
> flush(stderr()); flush(stdout())
>
> ### Name: Exponential
> ### Title: The Exponential Distribution
> ### Aliases: Exponential dexp pexp qexp rexp
> ### Keywords: distribution
>
> ### ** Examples
>
> dexp(1) - exp(-1) #-> 0
[1] 0>
>
>
> cleanEx(); ..nameEx <- "Fdist"
>
> ### * Fdist
>
> flush(stderr()); flush(stdout())
>
> ### Name: FDist
> ### Title: The F Distribution
> ### Aliases: FDist df pf qf rf
> ### Keywords: distribution
>
> ### ** Examples
>
> ## the density of the square of a t_m is 2*dt(x, m)/(2*x)
> # check this is the same as the density of F_{1,m}
> x <- seq(0.001, 5, len=100)
> all.equal(df(x^2, 1, 5), dt(x, 5)/x)
[1] TRUE>
> ## Identity: qf(2*p - 1, 1, df)) == qt(p, df)^2) for p >= 1/2
> p <- seq(1/2, .99, length=50); df <- 10
> rel.err <- function(x,y) ifelse(x==y,0, abs(x-y)/mean(abs(c(x,y))))
> quantile(rel.err(qf(2*p - 1, df1=1, df2=df), qt(p, df)^2), .90)# ~= 7e-9
90%
1.542479e-15 >
>
>
> cleanEx(); ..nameEx <- "GammaDist"
>
> ### * GammaDist
>
> flush(stderr()); flush(stdout())
>
> ### Name: GammaDist
> ### Title: The Gamma Distribution
> ### Aliases: GammaDist dgamma pgamma qgamma rgamma
> ### Keywords: distribution
>
> ### ** Examples
>
> -log(dgamma(1:4, shape=1))
[1] 1 2 3 4> p <- (1:9)/10
> pgamma(qgamma(p,shape=2), shape=2)
[1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9> 1 - 1/exp(qgamma(p, shape=1))
[1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9>
>
>
> cleanEx(); ..nameEx <- "Geometric"
>
> ### * Geometric
>
> flush(stderr()); flush(stdout())
>
> ### Name: Geometric
> ### Title: The Geometric Distribution
> ### Aliases: Geometric dgeom pgeom qgeom rgeom
> ### Keywords: distribution
>
> ### ** Examples
>
> qgeom((1:9)/10, prob = .2)
[1] 0 0 1 2 3 4 5 7 10> Ni <- rgeom(20, prob = 1/4); table(factor(Ni, 0:max(Ni)))
0 1 2 3 4 5 6 7 8 9 10 11
5 3 3 1 2 2 0 1 0 1 1 1 >
>
>
> cleanEx(); ..nameEx <- "Harman23.cor"
>
> ### * Harman23.cor
>
> flush(stderr()); flush(stdout())
>
> ### Name: Harman23.cor
> ### Title: Harman Example 2.3
> ### Aliases: Harman23.cor
> ### Keywords: datasets
>
> ### ** Examples
>
> data(Harman23.cor)
> (Harman23.FA <- factanal(factors = 1, covmat = Harman23.cor))
Call:
factanal(factors = 1, covmat = Harman23.cor)
Uniquenesses:
height arm.span forearm lower.leg weight
0.158 0.135 0.190 0.187 0.760
bitro.diameter chest.girth chest.width
0.829 0.877 0.801
Loadings:
Factor1
height 0.918
arm.span 0.930
forearm 0.900
lower.leg 0.902
weight 0.490
bitro.diameter 0.413
chest.girth 0.351
chest.width 0.446
Factor1
SS loadings 4.064
Proportion Var 0.508
Test of the hypothesis that 1 factor is sufficient.
The chi square statistic is 611.44 on 20 degrees of freedom.
The p-value is 1.12e-116 > for(factors in 2:4) print(update(Harman23.FA, factors = factors))
Error in La.svd(B) : infinite or missing values in x
Execution halted
If there is any other information I can give you, please write and I'll
send
it.
thanks;
Madeleine