R-intro.Rout.save

R version 2.11.0 Under development (unstable) (2010-01-07 r50925)
Copyright (C) 2010 The R Foundation for Statistical Computing
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 an HTML browser interface to help.
Type 'q()' to quit R.

> ####  "All the examples" from  ./R-intro.texi
> ####  -- in a way that this should be(come) an executable script.
> 
> options(digits=5, width=65)##--- for outputs !
> options(stringsAsFactors=TRUE) ## factory-fresh defaults
> options(useFancyQuotes=FALSE) ## avoid problems on Windows
> 
> ## 2. Simple Manipulations
> 
> x <- c(10.4, 5.6, 3.1, 6.4, 21.7)
> assign("x", c(10.4, 5.6, 3.1, 6.4, 21.7))
> c(10.4, 5.6, 3.1, 6.4, 21.7) -> x
> .Last.value
[1] 10.4  5.6  3.1  6.4 21.7
> 1/x
[1] 0.096154 0.178571 0.322581 0.156250 0.046083
> 
> y <- c(x, 0, x)
> v <- 2*x + y + 1
Warning message:
In 2 * x + y :
  longer object length is not a multiple of shorter object length
> ##-  Warning message:
> ##-  longer object length
> ##-  is not a multiple of shorter object length in: 2 * x + y
> 
> sqrt(-17)
[1] NaN
Warning message:
In sqrt(-17) : NaNs produced
> ##- [1] NaN
> ##- Warning message:
> ##- NaNs produced in: sqrt(-17)
> 
> sqrt(-17+0i)
[1] 0+4.1231i
> 
> ###--  2.3 .. regular sequences
> 
> 1:30
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
[21] 21 22 23 24 25 26 27 28 29 30
> 
> n <- 10
> 
> 1:n-1
 [1] 0 1 2 3 4 5 6 7 8 9
> 1:(n-1)
[1] 1 2 3 4 5 6 7 8 9
> 30:1
 [1] 30 29 28 27 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11
[21] 10  9  8  7  6  5  4  3  2  1
> 
> seq(2,10)
[1]  2  3  4  5  6  7  8  9 10
> all(seq(1,30) == seq(to=30, from=1))
[1] TRUE
> 
> seq(-5, 5, by=.2) -> s3
> s4 <- seq(length=51, from=-5, by=.2)
> all.equal(s3,s4)
[1] TRUE
> 
> s5 <- rep(x, times=5)
> s6 <- rep(x, each=5)
> 
> temp <- x > 13
> 
> z <- c(1:3,NA);  ind <- is.na(z)
> 
> 0/0
[1] NaN
> Inf - Inf
[1] NaN
> 
> labs <- paste(c("X","Y"), 1:10, sep="")
> labs
 [1] "X1"  "Y2"  "X3"  "Y4"  "X5"  "Y6"  "X7"  "Y8"  "X9"  "Y10"
> 
> x <- c(z,z-2)#-- NOT in texi ; more interesting
> y <- x[!is.na(x)]
> 
> (x+1)[(!is.na(x)) & x>0] -> z
> z
[1] 2 3 4 2
> 
> x <- c(x, 9:12)# long enough:
> x[1:10]
 [1]  1  2  3 NA -1  0  1 NA  9 10
> 
> c("x","y")[rep(c(1,2,2,1), times=4)]
 [1] "x" "y" "y" "x" "x" "y" "y" "x" "x" "y" "y" "x" "x" "y" "y"
[16] "x"
> 
> y <- x[-(1:5)]
> y
[1]  0  1 NA  9 10 11 12
> 
> fruit <- c(5, 10, 1, 20)
> names(fruit) <- c("orange", "banana", "apple", "peach")
>  fruit
orange banana  apple  peach 
     5     10      1     20 
> 
> lunch <- fruit[c("apple","orange")]
>  lunch
 apple orange 
     1      5 
> 
>  x
 [1]  1  2  3 NA -1  0  1 NA  9 10 11 12
> x[is.na(x)] <- 0
>  x
 [1]  1  2  3  0 -1  0  1  0  9 10 11 12
> 
> y <- -4:9
> y[y < 0] <- -y[y < 0]
> all(y == abs(y))
[1] TRUE
> y
 [1] 4 3 2 1 0 1 2 3 4 5 6 7 8 9
> 
> ###---------------
> 
> z <- 0:9
> digits <- as.character(z)
> digits
 [1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9"
> d <- as.integer(digits)
> all.equal(z, d)
[1] TRUE
> 
>  e <- numeric()
>  e[3] <- 17
>  e
[1] NA NA 17
> 
> alpha <- 10*(1:10)
> alpha <- alpha[2 * 1:5]
> alpha
[1]  20  40  60  80 100
> 
> winter <- data.frame(temp = c(-1,3,2,-2), cat = rep(c("A","B"), 2))
> winter
  temp cat
1   -1   A
2    3   B
3    2   A
4   -2   B
> unclass(winter)
$temp
[1] -1  3  2 -2

$cat
[1] A B A B
Levels: A B

attr(,"row.names")
[1] 1 2 3 4
> 
> ###------------ Ordered and unordered factors --------
> state <- c("tas", "sa",  "qld", "nsw", "nsw", "nt",  "wa",  "wa",
+            "qld", "vic", "nsw", "vic", "qld", "qld", "sa",  "tas",
+            "sa",  "nt",  "wa",  "vic", "qld", "nsw", "nsw", "wa",
+            "sa",  "act", "nsw", "vic", "vic", "act")
> statef <- factor(state)
> statef
 [1] tas sa  qld nsw nsw nt  wa  wa  qld vic nsw vic qld qld sa 
[16] tas sa  nt  wa  vic qld nsw nsw wa  sa  act nsw vic vic act
Levels: act nsw nt qld sa tas vic wa
> 
> levels(statef)
[1] "act" "nsw" "nt"  "qld" "sa"  "tas" "vic" "wa" 
> 
> incomes <- c(60, 49, 40, 61, 64, 60, 59, 54, 62, 69, 70, 42, 56,
+              61, 61, 61, 58, 51, 48, 65, 49, 49, 41, 48, 52, 46,
+              59, 46, 58, 43)
> 
> incmeans <- tapply(incomes, statef, mean)
> incmeans
   act    nsw     nt    qld     sa    tas    vic     wa 
44.500 57.333 55.500 53.600 55.000 60.500 56.000 52.250 
> 
> stderr <- function(x) sqrt(var(x)/length(x))
> 
> incster <- tapply(incomes, statef, stderr)
> incster
   act    nsw     nt    qld     sa    tas    vic     wa 
1.5000 4.3102 4.5000 4.1061 2.7386 0.5000 5.2440 2.6575 
> 
> ##
> z <- 1:1500
> dim(z) <- c(3,5,100)
> 
> 
> x <- array(1:20,dim=c(4,5))   # Generate a 4 by 5 array.
> x
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
> 
> i <- array(c(1:3,3:1),dim=c(3,2))
> i                             # @code{i} is a 3 by 2 index array.
     [,1] [,2]
[1,]    1    3
[2,]    2    2
[3,]    3    1
> 
> x[i]                          # Extract those elements
[1] 9 6 3
> 
> x[i] <- 0                     # Replace those elements by zeros.
> x
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    0   13   17
[2,]    2    0   10   14   18
[3,]    0    7   11   15   19
[4,]    4    8   12   16   20
> 
> n <- 60
> b <- 5 ; blocks    <- rep(1:b, length= n)
> v <- 6 ; varieties <- gl(v,10)
> 
> Xb <- matrix(0, n, b)
> Xv <- matrix(0, n, v)
> ib <- cbind(1:n, blocks)
> iv <- cbind(1:n, varieties)
> Xb[ib] <- 1
> Xv[iv] <- 1
> X <- cbind(Xb, Xv)
> 
> N <- crossprod(Xb, Xv)
> table(blocks,varieties)
      varieties
blocks 1 2 3 4 5 6
     1 2 2 2 2 2 2
     2 2 2 2 2 2 2
     3 2 2 2 2 2 2
     4 2 2 2 2 2 2
     5 2 2 2 2 2 2
> all(N == table(blocks,varieties))
[1] TRUE
> 
> h <- 1:17
> Z <- array(h, dim=c(3,4,2))
> ## If the size of  'h' is exactly 24
> h <- rep(h, length = 24)
> Z. <- Z ## the result is the same as
> Z <- h; dim(Z) <- c(3,4,2)
> stopifnot(identical(Z., Z))
> 
> Z <- array(0, c(3,4,2))
> 
> ## So if @code{A}, @code{B} and @code{C} are all similar arrays
> ## <init>
> A <- matrix(1:6, 3,2)
> B <- cbind(1, 1:3)
> C <- rbind(1, rbind(2, 3:4))
> stopifnot(dim(A) == dim(B),
+           dim(B) == dim(C))
> ## <init/>
> D <- 2*A*B + C + 1
> 
> a <- 1:9
> b <- 10*(1:3)
> 
> ab <- a %o% b
> stopifnot(ab == outer(a,b,"*"),
+           ab == outer(a,b))
> 
> x <- 1:10
> y <- -2:2
> f <- function(x, y) cos(y)/(1 + x^2)
> z <- outer(x, y, f)
> 
> 
> d <- outer(0:9, 0:9)
> fr <- table(outer(d, d, "-"))
> plot(as.numeric(names(fr)), fr, type="h",
+      xlab="Determinant", ylab="Frequency")
> 
> ##
> 
> B <- aperm(A, c(2,1))
> stopifnot(identical(B, t(A)))
> 
> ## for example, @code{A} and @code{B} are square matrices of the same size
> ## <init>
> A <- matrix(1:4, 2,2)
> B <- A - 1
> ## <init/>
> 
> A * B
     [,1] [,2]
[1,]    0    6
[2,]    2   12
> 
> A %*% B
     [,1] [,2]
[1,]    3   11
[2,]    4   16
> 
> 
> ## <init>
> x <- c(-1, 2)
> ## <init/>
> x %*% A %*% x
     [,1]
[1,]    7
> 
> x %*% x
     [,1]
[1,]    5
> stopifnot(x %*% x == sum(x^2))
> 
> xxT <- cbind(x) %*% x
> xxT
     [,1] [,2]
[1,]    1   -2
[2,]   -2    4
> stopifnot(identical(xxT, x %*% rbind(x)))
> 
> ## crossprod  ... (ADD)
> 
> ## diag  ... (ADD)
> 
> ## linear equations  ... (ADD)
> 
> ## solve  ... (ADD)
> 
> ## eigen:
> ## a symmetric matrix @code{Sm}
> ## <init>
> Sm <- matrix(-2:6, 3); Sm <- (Sm + t(Sm))/4; Sm
     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]    0    1    2
[3,]    1    2    3
> ## </init>
> ev <- eigen(Sm)
> 
> evals <- eigen(Sm)$values
> 
> ##  SVD .....
> 
> ## "if M is in fact square, then, ..."
> ## <init>
> M <- cbind(1,1:3,c(5,2,3))
> X <- cbind(1:9, .25*(-4:4)^2)
> X1 <- cbind(1:7, -1)
> X2 <- cbind(0,2:8)
> y <- c(1:4, 2:6)
> ## </init>
> 
> absdetM <- prod(svd(M)$d)
> stopifnot(all.equal(absdetM, abs(det(M))))# since det() nowadays exists
> 
> ans <- lsfit(X, y)
> 
>  Xplus <- qr(X)
>  b <- qr.coef(Xplus, y)
>  fit <- qr.fitted(Xplus, y)
>  res <- qr.resid(Xplus, y)
> ##
> 
> X <- cbind(1, X1, X2)
> 
> vec <- as.vector(X)
> vec <- c(X)
> 
> statefr <- table(statef)
> statefr
statef
act nsw  nt qld  sa tas vic  wa 
  2   6   2   5   4   2   5   4 
> statefr <- tapply(statef, statef, length)
> statefr
act nsw  nt qld  sa tas vic  wa 
  2   6   2   5   4   2   5   4 
> 
> factor(cut(incomes, breaks = 35+10*(0:7))) -> incomef
> table(incomef,statef)
         statef
incomef   act nsw nt qld sa tas vic wa
  (35,45]   1   1  0   1  0   0   1  0
  (45,55]   1   1  1   1  2   0   1  3
  (55,65]   0   3  1   3  2   2   2  1
  (65,75]   0   1  0   0  0   0   1  0
> 
> ###--- @chapter 6. Lists and data frames
> 
> Lst <- list(name="Fred", wife="Mary", no.children=3,
+             child.ages=c(4,7,9))
> Lst$name
[1] "Fred"
> Lst$wife
[1] "Mary"
> Lst$child.ages[1]
[1] 4
> stopifnot(Lst$name == Lst[[1]], Lst[[1]] == "Fred",
+           Lst$child.ages[1] == Lst[[4]][1], Lst[[4]][1] == 4
+           )
> 
> x <- "name" ; Lst[[x]]
[1] "Fred"
> 
> ## @section  6.2  Constructing and modifying lists
> 
> ##<init>
> Mat <- cbind(1, 2:4)
> ##</init>
> Lst[5] <- list(matrix=Mat)
> 
> ## @section  6.3  Data frames
> 
> accountants <- data.frame(home=statef, loot=incomes, shot=incomef)
> ## MM: add the next lines to R-intro.texi !
> accountants
   home loot    shot
1   tas   60 (55,65]
2    sa   49 (45,55]
3   qld   40 (35,45]
4   nsw   61 (55,65]
5   nsw   64 (55,65]
6    nt   60 (55,65]
7    wa   59 (55,65]
8    wa   54 (45,55]
9   qld   62 (55,65]
10  vic   69 (65,75]
11  nsw   70 (65,75]
12  vic   42 (35,45]
13  qld   56 (55,65]
14  qld   61 (55,65]
15   sa   61 (55,65]
16  tas   61 (55,65]
17   sa   58 (55,65]
18   nt   51 (45,55]
19   wa   48 (45,55]
20  vic   65 (55,65]
21  qld   49 (45,55]
22  nsw   49 (45,55]
23  nsw   41 (35,45]
24   wa   48 (45,55]
25   sa   52 (45,55]
26  act   46 (45,55]
27  nsw   59 (55,65]
28  vic   46 (45,55]
29  vic   58 (55,65]
30  act   43 (35,45]
> str(accountants)
'data.frame':	30 obs. of  3 variables:
 $ home: Factor w/ 8 levels "act","nsw","nt",..: 6 5 4 2 2 3 8 8 4 7 ...
 $ loot: num  60 49 40 61 64 60 59 54 62 69 ...
 $ shot: Factor w/ 4 levels "(35,45]","(45,55]",..: 3 2 1 3 3 3 3 2 3 4 ...
> 
> ## ..........
> 
> ###--- @chapter 8. Probability distributions
> 
> ## 2-tailed p-value for t distribution
> 2*pt(-2.43, df = 13)
[1] 0.030331
> ## upper 1% point for an F(2, 7)  distribution
> qf(0.01, 2, 7, lower.tail = FALSE)
[1] 9.5466
> 
> attach(faithful)
> summary(eruptions)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.60    2.16    4.00    3.49    4.45    5.10 
> 
> fivenum(eruptions)
[1] 1.6000 2.1585 4.0000 4.4585 5.1000
> 
> stem(eruptions)

  The decimal point is 1 digit(s) to the left of the |

  16 | 070355555588
  18 | 000022233333335577777777888822335777888
  20 | 00002223378800035778
  22 | 0002335578023578
  24 | 00228
  26 | 23
  28 | 080
  30 | 7
  32 | 2337
  34 | 250077
  36 | 0000823577
  38 | 2333335582225577
  40 | 0000003357788888002233555577778
  42 | 03335555778800233333555577778
  44 | 02222335557780000000023333357778888
  46 | 0000233357700000023578
  48 | 00000022335800333
  50 | 0370

> 
> hist(eruptions)
> 
> ## <IMG> postscript("images/hist.eps", ...)
> # make the bins smaller, make a plot of density
> hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)
> lines(density(eruptions, bw=0.1))
> rug(eruptions) # show the actual data points
> ## dev.off() <IMG/>
> 
> plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)
> 
> ## <IMG> postscript("images/ecdf.eps", ...)
> long <- eruptions[eruptions > 3]
> plot(ecdf(long), do.points=FALSE, verticals=TRUE)
> x <- seq(3, 5.4, 0.01)
> lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)
> ## dev.off() <IMG/>
> 
> par(pty="s")       # arrange for a square figure region
> qqnorm(long); qqline(long)
> 
> x <- rt(250, df = 5)
> qqnorm(x); qqline(x)
> 
> qqplot(qt(ppoints(250), df = 5), x, xlab = "Q-Q plot for t dsn")
> qqline(x)
> 
> shapiro.test(long)

	Shapiro-Wilk normality test

data:  long 
W = 0.9793, p-value = 0.01052

> 
> ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long)))

	One-sample Kolmogorov-Smirnov test

data:  long 
D = 0.0661, p-value = 0.4284
alternative hypothesis: two-sided 

Warning message:
In ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long))) :
  cannot compute correct p-values with ties
> 
> ##@section One- and two-sample tests
> 
> ## scan() from stdin :
> ## can be cut & pasted, but not parsed and hence not source()d
> ##scn A <- scan()
> ##scn 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97
> ##scn 80.05 80.03 80.02 80.00 80.02
> A <- c(79.98, 80.04, 80.02, 80.04, 80.03, 80.03, 80.04, 79.97,
+        80.05, 80.03, 80.02, 80, 80.02)
> ##scn B <- scan()
> ##scn 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97
> B <- c(80.02, 79.94, 79.98, 79.97, 79.97, 80.03, 79.95, 79.97)
> 
> ## <IMG> postscript("images/ice.eps", ...)
> boxplot(A, B)
> ## dev.off() <IMG/>
> 
> t.test(A, B)

	Welch Two Sample t-test

data:  A and B 
t = 3.2499, df = 12.027, p-value = 0.00694
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.013855 0.070183 
sample estimates:
mean of x mean of y 
   80.021    79.979 

> 
> var.test(A, B)

	F test to compare two variances

data:  A and B 
F = 0.5837, num df = 12, denom df = 7, p-value = 0.3938
alternative hypothesis: true ratio of variances is not equal to 1 
95 percent confidence interval:
 0.12511 2.10527 
sample estimates:
ratio of variances 
           0.58374 

> 
> t.test(A, B, var.equal=TRUE)

	Two Sample t-test

data:  A and B 
t = 3.4722, df = 19, p-value = 0.002551
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.016691 0.067348 
sample estimates:
mean of x mean of y 
   80.021    79.979 

> 
> wilcox.test(A, B)

	Wilcoxon rank sum test with continuity correction

data:  A and B 
W = 89, p-value = 0.007497
alternative hypothesis: true location shift is not equal to 0 

Warning message:
In wilcox.test.default(A, B) : cannot compute exact p-value with ties
> 
> plot(ecdf(A), do.points=FALSE, verticals=TRUE, xlim=range(A, B))
> plot(ecdf(B), do.points=FALSE, verticals=TRUE, add=TRUE)
> 
> ###--- @chapter Grouping, loops and conditional execution
> 
> 
> ###--- @chapter Writing your own functions
> 
> 
> ###--- @chapter Statistical models in R
> 
> 
> ###--- @chapter Graphical procedures
> 
> ###--- @appendix A sample session
> 
> ## "Simulate starting a new R session, by
> rm(list=ls(all=TRUE))
> set.seed(123) # for repeatability
> 
> if(interactive())
+     help.start()
> 
> x <- rnorm(50)
> y <- rnorm(x)
> plot(x, y)
> ls()
[1] "x" "y"
> rm(x, y)
> x <- 1:20
> w <- 1 + sqrt(x)/2
> dummy <- data.frame(x = x, y = x + rnorm(x)*w)
> dummy
    x        y
1   1 -0.06561
2   2  2.43853
3   3  2.53967
4   4  3.30491
5   5  2.98444
6   6  5.89982
7   7  5.17676
8   8  3.97323
9   9  8.04943
10 10 12.37206
11 11  9.47055
12 12 13.66099
13 13  8.46544
14 14 13.84049
15 15 16.52523
16 16 16.90346
17 17 17.32353
18 18 16.00015
19 19 16.29841
20 20 16.68585
> fm <- lm(y ~ x, data=dummy)
> summary(fm)

Call:
lm(formula = y ~ x, data = dummy)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5401 -1.1026 -0.0541  1.1525  3.2623 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.5431     0.8902   -0.61     0.55    
x             0.9653     0.0743   12.99  1.4e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.92 on 18 degrees of freedom
Multiple R-squared: 0.904,	Adjusted R-squared: 0.898 
F-statistic:  169 on 1 and 18 DF,  p-value: 1.39e-10 

> fm1 <- lm(y ~ x, data=dummy, weight=1/w^2)
> summary(fm1)

Call:
lm(formula = y ~ x, data = dummy, weights = 1/w^2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.32054 -0.44924 -0.00879  0.50876  1.26555 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.6155     0.6513   -0.94     0.36    
x             0.9721     0.0664   14.64  1.9e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.72 on 18 degrees of freedom
Multiple R-squared: 0.922,	Adjusted R-squared: 0.918 
F-statistic:  214 on 1 and 18 DF,  p-value: 1.94e-11 

> attach(dummy)
The following object(s) are masked _by_ '.GlobalEnv':

    x
> lrf <- lowess(x, y)
> plot(x, y)
> lines(x, lrf$y)
> abline(0, 1, lty=3)
> abline(coef(fm))
> abline(coef(fm1), col = "red")
> detach()# dummy
> 
> plot(fitted(fm), resid(fm),
+      xlab="Fitted values",
+      ylab="Residuals",
+      main="Residuals vs Fitted")
> qqnorm(resid(fm), main="Residuals Rankit Plot")
> rm(fm, fm1, lrf, x, dummy)
> 
> 
> filepath <- system.file("data", "morley.tab" , package="datasets")
> if(interactive()) file.show(filepath)
> mm <- read.table(filepath)
> mm
    Expt Run Speed
001    1   1   850
002    1   2   740
003    1   3   900
004    1   4  1070
005    1   5   930
006    1   6   850
007    1   7   950
008    1   8   980
009    1   9   980
010    1  10   880
011    1  11  1000
012    1  12   980
013    1  13   930
014    1  14   650
015    1  15   760
016    1  16   810
017    1  17  1000
018    1  18  1000
019    1  19   960
020    1  20   960
021    2   1   960
022    2   2   940
023    2   3   960
024    2   4   940
025    2   5   880
026    2   6   800
027    2   7   850
028    2   8   880
029    2   9   900
030    2  10   840
031    2  11   830
032    2  12   790
033    2  13   810
034    2  14   880
035    2  15   880
036    2  16   830
037    2  17   800
038    2  18   790
039    2  19   760
040    2  20   800
041    3   1   880
042    3   2   880
043    3   3   880
044    3   4   860
045    3   5   720
046    3   6   720
047    3   7   620
048    3   8   860
049    3   9   970
050    3  10   950
051    3  11   880
052    3  12   910
053    3  13   850
054    3  14   870
055    3  15   840
056    3  16   840
057    3  17   850
058    3  18   840
059    3  19   840
060    3  20   840
061    4   1   890
062    4   2   810
063    4   3   810
064    4   4   820
065    4   5   800
066    4   6   770
067    4   7   760
068    4   8   740
069    4   9   750
070    4  10   760
071    4  11   910
072    4  12   920
073    4  13   890
074    4  14   860
075    4  15   880
076    4  16   720
077    4  17   840
078    4  18   850
079    4  19   850
080    4  20   780
081    5   1   890
082    5   2   840
083    5   3   780
084    5   4   810
085    5   5   760
086    5   6   810
087    5   7   790
088    5   8   810
089    5   9   820
090    5  10   850
091    5  11   870
092    5  12   870
093    5  13   810
094    5  14   740
095    5  15   810
096    5  16   940
097    5  17   950
098    5  18   800
099    5  19   810
100    5  20   870
> mm$Expt <- factor(mm$Expt)
> mm$Run <- factor(mm$Run)
> attach(mm)
> plot(Expt, Speed, main="Speed of Light Data", xlab="Experiment No.")
> fm <- aov(Speed ~ Run + Expt, data=mm)
> summary(fm)
            Df Sum Sq Mean Sq F value Pr(>F)   
Run         19 113344    5965    1.11 0.3632   
Expt         4  94514   23629    4.38 0.0031 **
Residuals   76 410166    5397                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> fm0 <- update(fm, . ~ . - Run)
> anova(fm0, fm)
Analysis of Variance Table

Model 1: Speed ~ Expt
Model 2: Speed ~ Run + Expt
  Res.Df    RSS Df Sum of Sq    F Pr(>F)
1     95 523510                         
2     76 410166 19    113344 1.11   0.36
> detach()
> rm(fm, fm0)
> 
> x <- seq(-pi, pi, len=50)
> y <- x
> f <- outer(x, y, function(x, y) cos(y)/(1 + x^2))
> oldpar <- par(no.readonly = TRUE)
> par(pty="s")
> contour(x, y, f)
> contour(x, y, f, nlevels=15, add=TRUE)
> fa <- (f-t(f))/2
> contour(x, y, fa, nlevels=15)
> par(oldpar)
> image(x, y, f)
> image(x, y, fa)
> objects(); rm(x, y, f, fa)
[1] "f"        "fa"       "filepath" "mm"       "oldpar"  
[6] "w"        "x"        "y"       
> th <- seq(-pi, pi, len=100)
> z <- exp(1i*th)
> par(pty="s")
> plot(z, type="l")
> w <- rnorm(100) + rnorm(100)*1i
> w <- ifelse(Mod(w) > 1, 1/w, w)
> plot(w, xlim=c(-1,1), ylim=c(-1,1), pch="+",xlab="x", ylab="y")
> lines(z)
> 
> w <- sqrt(runif(100))*exp(2*pi*runif(100)*1i)
> plot(w, xlim=c(-1,1), ylim=c(-1,1), pch="+", xlab="x", ylab="y")
> lines(z)
> 
> rm(th, w, z)
> ## q()
> 
>