"Fossies" - the Fresh Open Source Software Archive

Member "R-3.6.0/doc/manual/R-intro.R" (25 Sep 2018, 9394 Bytes) of package /linux/misc/R-3.6.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) R source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 ####  "All the examples" from  ./R-intro.texi
    2 ####  -- in a way that this should be(come) an executable script.
    3 
    4 options(digits=5, width=65)##--- for outputs !
    5 options(stringsAsFactors=TRUE) ## factory-fresh defaults
    6 options(useFancyQuotes=FALSE) ## avoid problems on Windows
    7 
    8 ## 2. Simple Manipulations
    9 
   10 x <- c(10.4, 5.6, 3.1, 6.4, 21.7)
   11 assign("x", c(10.4, 5.6, 3.1, 6.4, 21.7))
   12 c(10.4, 5.6, 3.1, 6.4, 21.7) -> x
   13 .Last.value
   14 1/x
   15 
   16 y <- c(x, 0, x)
   17 v <- 2*x + y + 1
   18 ##-  Warning message:
   19 ##-  longer object length
   20 ##-  is not a multiple of shorter object length in: 2 * x + y
   21 
   22 sqrt(-17)
   23 ##- [1] NaN
   24 ##- Warning message:
   25 ##- NaNs produced in: sqrt(-17)
   26 
   27 sqrt(-17+0i)
   28 
   29 ###--  2.3 .. regular sequences
   30 
   31 1:30
   32 
   33 n <- 10
   34 
   35 1:n-1
   36 1:(n-1)
   37 30:1
   38 
   39 seq(2,10)
   40 all(seq(1,30) == seq(to=30, from=1))
   41 
   42 seq(-5, 5, by=.2) -> s3
   43 s4 <- seq(length=51, from=-5, by=.2)
   44 all.equal(s3,s4)
   45 
   46 s5 <- rep(x, times=5)
   47 s6 <- rep(x, each=5)
   48 
   49 temp <- x > 13
   50 
   51 z <- c(1:3,NA);  ind <- is.na(z)
   52 
   53 0/0
   54 Inf - Inf
   55 
   56 labs <- paste(c("X","Y"), 1:10, sep="")
   57 labs
   58 
   59 x <- c(z,z-2)#-- NOT in texi ; more interesting
   60 y <- x[!is.na(x)]
   61 
   62 (x+1)[(!is.na(x)) & x>0] -> z
   63 z
   64 
   65 x <- c(x, 9:12)# long enough:
   66 x[1:10]
   67 
   68 c("x","y")[rep(c(1,2,2,1), times=4)]
   69 
   70 y <- x[-(1:5)]
   71 y
   72 
   73 fruit <- c(5, 10, 1, 20)
   74 names(fruit) <- c("orange", "banana", "apple", "peach")
   75  fruit
   76 
   77 lunch <- fruit[c("apple","orange")]
   78  lunch
   79 
   80  x
   81 x[is.na(x)] <- 0
   82  x
   83 
   84 y <- -4:9
   85 y[y < 0] <- -y[y < 0]
   86 all(y == abs(y))
   87 y
   88 
   89 ###---------------
   90 
   91 z <- 0:9
   92 digits <- as.character(z)
   93 digits
   94 d <- as.integer(digits)
   95 all.equal(z, d)
   96 
   97  e <- numeric()
   98  e[3] <- 17
   99  e
  100 
  101 alpha <- 10*(1:10)
  102 alpha <- alpha[2 * 1:5]
  103 alpha
  104 
  105 winter <- data.frame(temp = c(-1,3,2,-2), cat = rep(c("A","B"), 2))
  106 winter
  107 unclass(winter)
  108 
  109 ###------------ Ordered and unordered factors --------
  110 state <- c("tas", "sa",  "qld", "nsw", "nsw", "nt",  "wa",  "wa",
  111            "qld", "vic", "nsw", "vic", "qld", "qld", "sa",  "tas",
  112            "sa",  "nt",  "wa",  "vic", "qld", "nsw", "nsw", "wa",
  113            "sa",  "act", "nsw", "vic", "vic", "act")
  114 statef <- factor(state)
  115 statef
  116 
  117 levels(statef)
  118 
  119 incomes <- c(60, 49, 40, 61, 64, 60, 59, 54, 62, 69, 70, 42, 56,
  120              61, 61, 61, 58, 51, 48, 65, 49, 49, 41, 48, 52, 46,
  121              59, 46, 58, 43)
  122 
  123 incmeans <- tapply(incomes, statef, mean)
  124 incmeans
  125 
  126 stdError <- function(x) sqrt(var(x)/length(x))
  127 
  128 incster <- tapply(incomes, statef, stdError)
  129 incster
  130 
  131 ##
  132 z <- 1:1500
  133 dim(z) <- c(3,5,100)
  134 
  135 
  136 x <- array(1:20,dim=c(4,5))   # Generate a 4 by 5 array.
  137 x
  138 
  139 i <- array(c(1:3,3:1),dim=c(3,2))
  140 i                             # @code{i} is a 3 by 2 index array.
  141 
  142 x[i]                          # Extract those elements
  143 
  144 x[i] <- 0                     # Replace those elements by zeros.
  145 x
  146 
  147 n <- 60
  148 b <- 5 ; blocks    <- rep(1:b, length= n)
  149 v <- 6 ; varieties <- gl(v,10)
  150 
  151 Xb <- matrix(0, n, b)
  152 Xv <- matrix(0, n, v)
  153 ib <- cbind(1:n, blocks)
  154 iv <- cbind(1:n, varieties)
  155 Xb[ib] <- 1
  156 Xv[iv] <- 1
  157 X <- cbind(Xb, Xv)
  158 
  159 N <- crossprod(Xb, Xv)
  160 table(blocks,varieties)
  161 all(N == table(blocks,varieties))
  162 
  163 h <- 1:17
  164 Z <- array(h, dim=c(3,4,2))
  165 ## If the size of  'h' is exactly 24
  166 h <- rep(h, length = 24)
  167 Z. <- Z ## the result is the same as
  168 Z <- h; dim(Z) <- c(3,4,2)
  169 stopifnot(identical(Z., Z))
  170 
  171 Z <- array(0, c(3,4,2))
  172 
  173 ## So if @code{A}, @code{B} and @code{C} are all similar arrays
  174 ## <init>
  175 A <- matrix(1:6, 3,2)
  176 B <- cbind(1, 1:3)
  177 C <- rbind(1, rbind(2, 3:4))
  178 stopifnot(dim(A) == dim(B),
  179           dim(B) == dim(C))
  180 ## <init/>
  181 D <- 2*A*B + C + 1
  182 
  183 a <- 1:9
  184 b <- 10*(1:3)
  185 
  186 ab <- a %o% b
  187 stopifnot(ab == outer(a,b,"*"),
  188           ab == outer(a,b))
  189 
  190 x <- 1:10
  191 y <- -2:2
  192 f <- function(x, y) cos(y)/(1 + x^2)
  193 z <- outer(x, y, f)
  194 
  195 
  196 d <- outer(0:9, 0:9)
  197 fr <- table(outer(d, d, "-"))
  198 plot(fr, xlab="Determinant", ylab="Frequency")
  199 
  200 ##
  201 
  202 B <- aperm(A, c(2,1))
  203 stopifnot(identical(B, t(A)))
  204 
  205 ## for example, @code{A} and @code{B} are square matrices of the same size
  206 ## <init>
  207 A <- matrix(1:4, 2,2)
  208 B <- A - 1
  209 ## <init/>
  210 
  211 A * B
  212 
  213 A %*% B
  214 
  215 
  216 ## <init>
  217 x <- c(-1, 2)
  218 ## <init/>
  219 x %*% A %*% x
  220 
  221 x %*% x
  222 stopifnot(x %*% x == sum(x^2))
  223 
  224 xxT <- cbind(x) %*% x
  225 xxT
  226 stopifnot(identical(xxT, x %*% rbind(x)))
  227 
  228 ## crossprod  ... (ADD)
  229 
  230 ## diag  ... (ADD)
  231 
  232 ## linear equations  ... (ADD)
  233 
  234 ## solve  ... (ADD)
  235 
  236 ## eigen:
  237 ## a symmetric matrix @code{Sm}
  238 ## <init>
  239 Sm <- matrix(-2:6, 3); Sm <- (Sm + t(Sm))/4; Sm
  240 ## </init>
  241 ev <- eigen(Sm)
  242 
  243 evals <- eigen(Sm)$values
  244 
  245 ##  SVD .....
  246 
  247 ## "if M is in fact square, then, ..."
  248 ## <init>
  249 M <- cbind(1,1:3,c(5,2,3))
  250 X <- cbind(1:9, .25*(-4:4)^2)
  251 X1 <- cbind(1:7, -1)
  252 X2 <- cbind(0,2:8)
  253 y <- c(1:4, 2:6)
  254 ## </init>
  255 
  256 absdetM <- prod(svd(M)$d)
  257 stopifnot(all.equal(absdetM, abs(det(M))))# since det() nowadays exists
  258 
  259 ans <- lsfit(X, y)
  260 
  261  Xplus <- qr(X)
  262  b <- qr.coef(Xplus, y)
  263  fit <- qr.fitted(Xplus, y)
  264  res <- qr.resid(Xplus, y)
  265 ##
  266 
  267 X <- cbind(1, X1, X2)
  268 
  269 vec <- as.vector(X)
  270 vec <- c(X)
  271 
  272 statefr <- table(statef)
  273 statefr
  274 statefr <- tapply(statef, statef, length)
  275 statefr
  276 
  277 factor(cut(incomes, breaks = 35+10*(0:7))) -> incomef
  278 table(incomef,statef)
  279 
  280 ###--- @chapter 6. Lists and data frames
  281 
  282 Lst <- list(name="Fred", wife="Mary", no.children=3,
  283             child.ages=c(4,7,9))
  284 Lst$name
  285 Lst$wife
  286 Lst$child.ages[1]
  287 stopifnot(Lst$name == Lst[[1]], Lst[[1]] == "Fred",
  288           Lst$child.ages[1] == Lst[[4]][1], Lst[[4]][1] == 4
  289           )
  290 
  291 x <- "name" ; Lst[[x]]
  292 
  293 ## @section  6.2  Constructing and modifying lists
  294 
  295 ##<init>
  296 Mat <- cbind(1, 2:4)
  297 ##</init>
  298 Lst[5] <- list(matrix=Mat)
  299 
  300 ## @section  6.3  Data frames
  301 
  302 accountants <- data.frame(home=statef, loot=incomes, shot=incomef)
  303 ## MM: add the next lines to R-intro.texi !
  304 accountants
  305 str(accountants)
  306 
  307 ## ..........
  308 
  309 ###--- @chapter 8. Probability distributions
  310 
  311 ## 2-tailed p-value for t distribution
  312 2*pt(-2.43, df = 13)
  313 ## upper 1% point for an F(2, 7)  distribution
  314 qf(0.01, 2, 7, lower.tail = FALSE)
  315 
  316 attach(faithful)
  317 summary(eruptions)
  318 
  319 fivenum(eruptions)
  320 
  321 stem(eruptions)
  322 
  323 hist(eruptions)
  324 
  325 ## <IMG> postscript("images/hist.eps", ...)
  326 # make the bins smaller, make a plot of density
  327 hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)
  328 lines(density(eruptions, bw=0.1))
  329 rug(eruptions) # show the actual data points
  330 ## dev.off() <IMG/>
  331 
  332 plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)
  333 
  334 ## <IMG> postscript("images/ecdf.eps", ...)
  335 long <- eruptions[eruptions > 3]
  336 plot(ecdf(long), do.points=FALSE, verticals=TRUE)
  337 x <- seq(3, 5.4, 0.01)
  338 lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)
  339 ## dev.off() <IMG/>
  340 
  341 par(pty="s")       # arrange for a square figure region
  342 qqnorm(long); qqline(long)
  343 
  344 x <- rt(250, df = 5)
  345 qqnorm(x); qqline(x)
  346 
  347 qqplot(qt(ppoints(250), df = 5), x, xlab = "Q-Q plot for t dsn")
  348 qqline(x)
  349 
  350 shapiro.test(long)
  351 
  352 ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long)))
  353 
  354 ##@section One- and two-sample tests
  355 
  356 ## scan() from stdin :
  357 ## can be cut & pasted, but not parsed and hence not source()d
  358 ##scn A <- scan()
  359 ##scn 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97
  360 ##scn 80.05 80.03 80.02 80.00 80.02
  361 A <- c(79.98, 80.04, 80.02, 80.04, 80.03, 80.03, 80.04, 79.97,
  362        80.05, 80.03, 80.02, 80, 80.02)
  363 ##scn B <- scan()
  364 ##scn 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97
  365 B <- c(80.02, 79.94, 79.98, 79.97, 79.97, 80.03, 79.95, 79.97)
  366 
  367 ## <IMG> postscript("images/ice.eps", ...)
  368 boxplot(A, B)
  369 ## dev.off() <IMG/>
  370 
  371 t.test(A, B)
  372 
  373 var.test(A, B)
  374 
  375 t.test(A, B, var.equal=TRUE)
  376 
  377 wilcox.test(A, B)
  378 
  379 plot(ecdf(A), do.points=FALSE, verticals=TRUE, xlim=range(A, B))
  380 plot(ecdf(B), do.points=FALSE, verticals=TRUE, add=TRUE)
  381 
  382 ###--- @chapter Grouping, loops and conditional execution
  383 
  384 
  385 ###--- @chapter Writing your own functions
  386 
  387 
  388 ###--- @chapter Statistical models in R
  389 
  390 
  391 ###--- @chapter Graphical procedures
  392 
  393 ###--- @appendix A sample session
  394 
  395 ## "Simulate starting a new R session, by
  396 rm(list=ls(all=TRUE))
  397 set.seed(123) # for repeatability
  398 
  399 if(interactive())
  400     help.start()
  401 
  402 x <- rnorm(50)
  403 y <- rnorm(x)
  404 plot(x, y)
  405 ls()
  406 rm(x, y)
  407 x <- 1:20
  408 w <- 1 + sqrt(x)/2
  409 dummy <- data.frame(x = x, y = x + rnorm(x)*w)
  410 dummy
  411 fm <- lm(y ~ x, data=dummy)
  412 summary(fm)
  413 fm1 <- lm(y ~ x, data=dummy, weight=1/w^2)
  414 summary(fm1)
  415 attach(dummy)
  416 lrf <- lowess(x, y)
  417 plot(x, y)
  418 lines(x, lrf$y)
  419 abline(0, 1, lty=3)
  420 abline(coef(fm))
  421 abline(coef(fm1), col = "red")
  422 detach()# dummy
  423 
  424 plot(fitted(fm), resid(fm),
  425      xlab="Fitted values",
  426      ylab="Residuals",
  427      main="Residuals vs Fitted")
  428 qqnorm(resid(fm), main="Residuals Rankit Plot")
  429 rm(fm, fm1, lrf, x, dummy)
  430 
  431 
  432 filepath <- system.file("data", "morley.tab" , package="datasets")
  433 if(interactive()) file.show(filepath)
  434 mm <- read.table(filepath)
  435 mm
  436 mm$Expt <- factor(mm$Expt)
  437 mm$Run <- factor(mm$Run)
  438 attach(mm)
  439 plot(Expt, Speed, main="Speed of Light Data", xlab="Experiment No.")
  440 fm <- aov(Speed ~ Run + Expt, data=mm)
  441 summary(fm)
  442 fm0 <- update(fm, . ~ . - Run)
  443 anova(fm0, fm)
  444 detach()
  445 rm(fm, fm0)
  446 
  447 x <- seq(-pi, pi, len=50)
  448 y <- x
  449 f <- outer(x, y, function(x, y) cos(y)/(1 + x^2))
  450 oldpar <- par(no.readonly = TRUE)
  451 par(pty="s")
  452 contour(x, y, f)
  453 contour(x, y, f, nlevels=15, add=TRUE)
  454 fa <- (f-t(f))/2
  455 contour(x, y, fa, nlevels=15)
  456 par(oldpar)
  457 image(x, y, f)
  458 image(x, y, fa)
  459 objects(); rm(x, y, f, fa)
  460 th <- seq(-pi, pi, len=100)
  461 z <- exp(1i*th)
  462 par(pty="s")
  463 plot(z, type="l")
  464 w <- rnorm(100) + rnorm(100)*1i
  465 w <- ifelse(Mod(w) > 1, 1/w, w)
  466 plot(w, xlim=c(-1,1), ylim=c(-1,1), pch="+",xlab="x", ylab="y")
  467 lines(z)
  468 
  469 w <- sqrt(runif(100))*exp(2*pi*runif(100)*1i)
  470 plot(w, xlim=c(-1,1), ylim=c(-1,1), pch="+", xlab="x", ylab="y")
  471 lines(z)
  472 
  473 rm(th, w, z)
  474 ## q()
  475