1 2R version 4.1.0 Patched (2021-06-17 r80511) -- "Camp Pontanezen" 3Copyright (C) 2021 The R Foundation for Statistical Computing 4Platform: x86_64-pc-linux-gnu (64-bit) 5 6R is free software and comes with ABSOLUTELY NO WARRANTY. 7You are welcome to redistribute it under certain conditions. 8Type 'license()' or 'licence()' for distribution details. 9 10R is a collaborative project with many contributors. 11Type 'contributors()' for more information and 12'citation()' on how to cite R or R packages in publications. 13 14Type 'demo()' for some demos, 'help()' for on-line help, or 15'help.start()' for an HTML browser interface to help. 16Type 'q()' to quit R. 17 18> #### Run all demos that do not depend on tcl and other specials. 19> .ptime <- proc.time() 20> set.seed(123) 21> options(keep.source=TRUE, useFancyQuotes=FALSE, warn = 1) 22> 23> ## Drop these for strict testing {and add them to demos2.R) 24> ## lm.glm is in ../src/library/utils/man/demo.Rd }: 25> dont <- list(graphics = c("Hershey", "Japanese", "plotmath"), 26+ stats = c("lm.glm", "nlm") 27+ ) 28> ## don't take tcltk here 29> for(pkg in c("base", "graphics", "stats")) { 30+ 31+ demos <- list.files(file.path(system.file(package = pkg), "demo"), 32+ pattern = "\\.R$") 33+ demos <- demos[is.na(match(demos, paste(dont[[pkg]], "R",sep=".")))] 34+ 35+ if(length(demos)) { 36+ if(need <- pkg != "base" && 37+ !any((fpkg <- paste("package", pkg, sep=":")) == search())) 38+ library(pkg, character.only = TRUE) 39+ 40+ for(nam in sub("\\.R$", "", demos)) 41+ demo(nam, character.only = TRUE) 42+ 43+ if(need) detach(pos = which(fpkg == search())) 44+ } 45+ } 46 47 48 demo(error.catching) 49 ---- ~~~~~~~~~~~~~~ 50 51> ##================================================================## 52> ### In longer simulations, aka computer experiments, ### 53> ### you may want to ### 54> ### 1) catch all errors and warnings (and continue) ### 55> ### 2) store the error or warning messages ### 56> ### ### 57> ### Here's a solution (see R-help mailing list, Dec 9, 2010): ### 58> ##================================================================## 59> 60> ##' Catch *and* save both errors and warnings, and in the case of 61> ##' a warning, also keep the computed result. 62> ##' 63> ##' @title tryCatch both warnings (with value) and errors 64> ##' @param expr an \R expression to evaluate 65> ##' @return a list with 'value' and 'warning', where 66> ##' 'value' may be an error caught. 67> ##' @author Martin Maechler; 68> ##' Copyright (C) 2010-2012 The R Core Team 69> tryCatch.W.E <- function(expr) 70+ { 71+ W <- NULL 72+ w.handler <- function(w){ # warning handler 73+ W <<- w 74+ invokeRestart("muffleWarning") 75+ } 76+ list(value = withCallingHandlers(tryCatch(expr, error = function(e) e), 77+ warning = w.handler), 78+ warning = W) 79+ } 80 81> str( tryCatch.W.E( log( 2 ) ) ) 82List of 2 83 $ value : num 0.693 84 $ warning: NULL 85 86> str( tryCatch.W.E( log( -1) ) ) 87List of 2 88 $ value : num NaN 89 $ warning:List of 2 90 ..$ message: chr "NaNs produced" 91 ..$ call : language log(-1) 92 ..- attr(*, "class")= chr [1:3] "simpleWarning" "warning" "condition" 93 94> str( tryCatch.W.E( log("a") ) ) 95List of 2 96 $ value :List of 2 97 ..$ message: chr "non-numeric argument to mathematical function" 98 ..$ call : language log("a") 99 ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition" 100 $ warning: NULL 101 102 103 demo(is.things) 104 ---- ~~~~~~~~~ 105 106> # Copyright (C) 1997-2018 The R Core Team 107> 108> ### The Base package has a couple of non-functions: 109> ## 110> ## These may be in "base" when they exist; discount them here 111> ## (see also 'dont.mind' in checkConflicts() inside library()) : 112> xtraBaseNms <- c("last.dump", "last.warning", ".Last.value", 113+ ".Random.seed", ".Traceback") 114 115> ls.base <- Filter(function(nm) is.na(match(nm, xtraBaseNms)), 116+ ls("package:base", all=TRUE)) 117 118> base.is.f <- sapply(ls.base, function(x) is.function(get(x))) 119 120> cat("\nNumber of all base objects:\t", length(ls.base), 121+ "\nNumber of functions from these:\t", sum(base.is.f), 122+ "\n\t starting with 'is.' :\t ", 123+ sum(grepl("^is\\.", ls.base[base.is.f])), "\n", sep = "") 124 125Number of all base objects: 1369 126Number of functions from these: 1324 127 starting with 'is.' : 50 128 129> ## R ver.| #{is*()} 130> ## ------+--------- 131> ## 0.14 : 31 132> ## 0.50 : 33 133> ## 0.60 : 34 134> ## 0.63 : 37 135> ## 1.0.0 : 38 136> ## 1.3.0 : 41 137> ## 1.6.0 : 45 138> ## 2.0.0 : 45 139> ## 2.7.0 : 48 140> ## 3.0.0 : 49 141> if(interactive()) { 142+ nonDots <- function(nm) substr(nm, 1L, 1L) != "." 143+ cat("Base non-functions not starting with \".\":\n") 144+ Filter(nonDots, ls.base[!base.is.f]) 145+ } 146 147> ## Do we have a method (probably)? 148> is.method <- function(fname) { 149+ isFun <- function(name) (exists(name, mode="function") && 150+ is.na(match(name, c("is", "as")))) 151+ np <- length(sp <- strsplit(fname, split = "\\.")[[1]]) 152+ if(np <= 1 ) 153+ FALSE 154+ else 155+ (isFun(paste(sp[1:(np-1)], collapse = '.')) || 156+ (np >= 3 && 157+ isFun(paste(sp[1:(np-2)], collapse = '.')))) 158+ } 159 160> is.ALL <- function(obj, func.names = ls(pos=length(search())), 161+ not.using = c("is.single", "is.real", "is.loaded", 162+ "is.empty.model", "is.R", "is.element", "is.unsorted"), 163+ true.only = FALSE, debug = FALSE) 164+ { 165+ ## Purpose: show many 'attributes' of R object __obj__ 166+ ## ------------------------------------------------------------------------- 167+ ## Arguments: obj: any R object 168+ ## ------------------------------------------------------------------------- 169+ ## Author: Martin Maechler, Date: 6 Dec 1996 170+ 171+ is.fn <- func.names[substring(func.names,1,3) == "is."] 172+ is.fn <- is.fn[substring(is.fn,1,7) != "is.na<-"] 173+ use.fn <- is.fn[ is.na(match(is.fn, not.using)) 174+ & ! sapply(is.fn, is.method) ] 175+ 176+ r <- if(true.only) character(0) 177+ else structure(vector("list", length= length(use.fn)), names= use.fn) 178+ for(f in use.fn) { 179+ if(any(f == c("is.na", "is.finite"))) { 180+ if(!is.list(obj) && !is.vector(obj) && !is.array(obj)) { 181+ if(!true.only) r[[f]] <- NA 182+ next 183+ } 184+ } 185+ if(any(f == c("is.nan", "is.finite", "is.infinite"))) { 186+ if(!is.atomic(obj)) { 187+ if(!true.only) r[[f]] <- NA 188+ next 189+ } 190+ } 191+ if(debug) cat(f,"") 192+ fn <- get(f) 193+ rr <- if(is.primitive(fn) || length(formals(fn))>0) fn(obj) else fn() 194+ if(!is.logical(rr)) cat("f=",f," --- rr is NOT logical = ",rr,"\n") 195+ ##if(1!=length(rr)) cat("f=",f," --- rr NOT of length 1; = ",rr,"\n") 196+ if(true.only && length(rr)==1 && !is.na(rr) && rr) r <- c(r, f) 197+ else if(!true.only) r[[f]] <- rr 198+ } 199+ if(debug)cat("\n") 200+ if(is.list(r)) structure(r, class = "isList") else r 201+ } 202 203> print.isList <- function(x, ..., verbose = getOption("verbose")) 204+ { 205+ ## Purpose: print METHOD for `isList' objects 206+ ## ------------------------------------------------ 207+ ## Author: Martin Maechler, Date: 12 Mar 1997 208+ if(is.list(x)) { 209+ if(verbose) cat("print.isList(): list case (length=",length(x),")\n") 210+ nm <- format(names(x)) 211+ rr <- lapply(x, stats::symnum, na = "NA") 212+ for(i in seq_along(x)) cat(nm[i],":",rr[[i]],"\n", ...) 213+ } else NextMethod("print", ...) 214+ } 215 216> is.ALL(NULL) 217is.array : . 218is.atomic : | 219is.call : . 220is.character : . 221is.complex : . 222is.data.frame : . 223is.double : . 224is.environment : . 225is.expression : . 226is.factor : . 227is.finite : NA 228is.function : . 229is.infinite : 230is.integer : . 231is.language : . 232is.list : . 233is.logical : . 234is.matrix : . 235is.na : NA 236is.name : . 237is.nan : 238is.null : | 239is.numeric : . 240is.numeric_version : . 241is.object : . 242is.ordered : . 243is.package_version : . 244is.pairlist : | 245is.primitive : . 246is.qr : . 247is.raw : . 248is.recursive : . 249is.symbol : . 250is.table : . 251is.vector : . 252 253> ##fails: is.ALL(NULL, not.using = c("is.single", "is.loaded")) 254> is.ALL(NULL, true.only = TRUE) 255[1] "is.atomic" "is.null" "is.pairlist" 256 257> all.equal(NULL, pairlist()) 258[1] TRUE 259 260> ## list() != NULL == pairlist() : 261> is.ALL(list(), true.only = TRUE) 262[1] "is.list" "is.recursive" "is.vector" 263 264> (pl <- is.ALL(pairlist(1, list(3,"A")), true.only = TRUE)) 265[1] "is.list" "is.pairlist" "is.recursive" 266 267> (ll <- is.ALL( list(1,pairlist(3,"A")), true.only = TRUE)) 268[1] "is.list" "is.recursive" "is.vector" 269 270> all.equal(pl[pl != "is.pairlist"], 271+ ll[ll != "is.vector"])## TRUE 272[1] TRUE 273 274> is.ALL(1:5) 275is.array : . 276is.atomic : | 277is.call : . 278is.character : . 279is.complex : . 280is.data.frame : . 281is.double : . 282is.environment : . 283is.expression : . 284is.factor : . 285is.finite : | | | | | 286is.function : . 287is.infinite : . . . . . 288is.integer : | 289is.language : . 290is.list : . 291is.logical : . 292is.matrix : . 293is.na : . . . . . 294is.name : . 295is.nan : . . . . . 296is.null : . 297is.numeric : | 298is.numeric_version : . 299is.object : . 300is.ordered : . 301is.package_version : . 302is.pairlist : . 303is.primitive : . 304is.qr : . 305is.raw : . 306is.recursive : . 307is.symbol : . 308is.table : . 309is.vector : | 310 311> is.ALL(array(1:24, 2:4)) 312is.array : | 313is.atomic : | 314is.call : . 315is.character : . 316is.complex : . 317is.data.frame : . 318is.double : . 319is.environment : . 320is.expression : . 321is.factor : . 322is.finite : | | | | | | | | | | | | | | | | | | | | | | | | 323is.function : . 324is.infinite : . . . . . . . . . . . . . . . . . . . . . . . . 325is.integer : | 326is.language : . 327is.list : . 328is.logical : . 329is.matrix : . 330is.na : . . . . . . . . . . . . . . . . . . . . . . . . 331is.name : . 332is.nan : . . . . . . . . . . . . . . . . . . . . . . . . 333is.null : . 334is.numeric : | 335is.numeric_version : . 336is.object : . 337is.ordered : . 338is.package_version : . 339is.pairlist : . 340is.primitive : . 341is.qr : . 342is.raw : . 343is.recursive : . 344is.symbol : . 345is.table : . 346is.vector : . 347 348> is.ALL(1 + 3) 349is.array : . 350is.atomic : | 351is.call : . 352is.character : . 353is.complex : . 354is.data.frame : . 355is.double : | 356is.environment : . 357is.expression : . 358is.factor : . 359is.finite : | 360is.function : . 361is.infinite : . 362is.integer : . 363is.language : . 364is.list : . 365is.logical : . 366is.matrix : . 367is.na : . 368is.name : . 369is.nan : . 370is.null : . 371is.numeric : | 372is.numeric_version : . 373is.object : . 374is.ordered : . 375is.package_version : . 376is.pairlist : . 377is.primitive : . 378is.qr : . 379is.raw : . 380is.recursive : . 381is.symbol : . 382is.table : . 383is.vector : | 384 385> e13 <- expression(1 + 3) 386 387> is.ALL(e13) 388Warning in fn(obj) : 389 is.na() applied to non-(list or vector) of type 'expression' 390is.array : . 391is.atomic : . 392is.call : . 393is.character : . 394is.complex : . 395is.data.frame : . 396is.double : . 397is.environment : . 398is.expression : | 399is.factor : . 400is.finite : NA 401is.function : . 402is.infinite : NA 403is.integer : . 404is.language : | 405is.list : . 406is.logical : . 407is.matrix : . 408is.na : . 409is.name : . 410is.nan : NA 411is.null : . 412is.numeric : . 413is.numeric_version : . 414is.object : . 415is.ordered : . 416is.package_version : . 417is.pairlist : . 418is.primitive : . 419is.qr : . 420is.raw : . 421is.recursive : | 422is.symbol : . 423is.table : . 424is.vector : | 425 426> is.ALL(substitute(expression(a + 3), list(a=1)), true.only = TRUE) 427[1] "is.call" "is.language" "is.recursive" 428 429> is.ALL(y ~ x) #--> NA for is.na & is.finite 430is.array : . 431is.atomic : . 432is.call : | 433is.character : . 434is.complex : . 435is.data.frame : . 436is.double : . 437is.environment : . 438is.expression : . 439is.factor : . 440is.finite : NA 441is.function : . 442is.infinite : NA 443is.integer : . 444is.language : | 445is.list : . 446is.logical : . 447is.matrix : . 448is.na : NA 449is.name : . 450is.nan : NA 451is.null : . 452is.numeric : . 453is.numeric_version : . 454is.object : | 455is.ordered : . 456is.package_version : . 457is.pairlist : . 458is.primitive : . 459is.qr : . 460is.raw : . 461is.recursive : | 462is.symbol : . 463is.table : . 464is.vector : . 465 466> is0 <- is.ALL(numeric(0)) 467 468> is0.ok <- 1 == (lis0 <- sapply(is0, length)) 469 470> is0[!is0.ok] 471$is.finite 472logical(0) 473 474$is.infinite 475logical(0) 476 477$is.na 478logical(0) 479 480$is.nan 481logical(0) 482 483 484> is0 <- unlist(is0) 485 486> is0 487 is.array is.atomic is.call is.character 488 FALSE TRUE FALSE FALSE 489 is.complex is.data.frame is.double is.environment 490 FALSE FALSE TRUE FALSE 491 is.expression is.factor is.function is.integer 492 FALSE FALSE FALSE FALSE 493 is.language is.list is.logical is.matrix 494 FALSE FALSE FALSE FALSE 495 is.name is.null is.numeric is.numeric_version 496 FALSE FALSE TRUE FALSE 497 is.object is.ordered is.package_version is.pairlist 498 FALSE FALSE FALSE FALSE 499 is.primitive is.qr is.raw is.recursive 500 FALSE FALSE FALSE FALSE 501 is.symbol is.table is.vector 502 FALSE FALSE TRUE 503 504> ispi <- unlist(is.ALL(pi)) 505 506> all(ispi[is0.ok] == is0) 507[1] TRUE 508 509> is.ALL(numeric(0), true=TRUE) 510[1] "is.atomic" "is.double" "is.numeric" "is.vector" 511 512> is.ALL(array(1,1:3), true=TRUE) 513[1] "is.array" "is.atomic" "is.double" "is.numeric" 514 515> is.ALL(cbind(1:3), true=TRUE) 516[1] "is.array" "is.atomic" "is.integer" "is.matrix" "is.numeric" 517 518> is.ALL(structure(1:7, names = paste("a",1:7,sep=""))) 519is.array : . 520is.atomic : | 521is.call : . 522is.character : . 523is.complex : . 524is.data.frame : . 525is.double : . 526is.environment : . 527is.expression : . 528is.factor : . 529is.finite : | | | | | | | 530is.function : . 531is.infinite : . . . . . . . 532is.integer : | 533is.language : . 534is.list : . 535is.logical : . 536is.matrix : . 537is.na : . . . . . . . 538is.name : . 539is.nan : . . . . . . . 540is.null : . 541is.numeric : | 542is.numeric_version : . 543is.object : . 544is.ordered : . 545is.package_version : . 546is.pairlist : . 547is.primitive : . 548is.qr : . 549is.raw : . 550is.recursive : . 551is.symbol : . 552is.table : . 553is.vector : | 554 555> is.ALL(structure(1:7, names = paste("a",1:7,sep="")), true.only = TRUE) 556[1] "is.atomic" "is.integer" "is.numeric" "is.vector" 557 558> x <- 1:20 ; y <- 5 + 6*x + rnorm(20) 559 560> lm.xy <- lm(y ~ x) 561 562> is.ALL(lm.xy) 563is.array : . 564is.atomic : . 565is.call : . 566is.character : . 567is.complex : . 568is.data.frame : . 569is.double : . 570is.environment : . 571is.expression : . 572is.factor : . 573is.finite : NA 574is.function : . 575is.infinite : NA 576is.integer : . 577is.language : . 578is.list : | 579is.logical : . 580is.matrix : . 581is.na : . . . . . . . . . . . . 582is.name : . 583is.nan : NA 584is.null : . 585is.numeric : . 586is.numeric_version : . 587is.object : | 588is.ordered : . 589is.package_version : . 590is.pairlist : . 591is.primitive : . 592is.qr : . 593is.raw : . 594is.recursive : | 595is.symbol : . 596is.table : . 597is.vector : . 598 599> is.ALL(structure(1:7, names = paste("a",1:7,sep=""))) 600is.array : . 601is.atomic : | 602is.call : . 603is.character : . 604is.complex : . 605is.data.frame : . 606is.double : . 607is.environment : . 608is.expression : . 609is.factor : . 610is.finite : | | | | | | | 611is.function : . 612is.infinite : . . . . . . . 613is.integer : | 614is.language : . 615is.list : . 616is.logical : . 617is.matrix : . 618is.na : . . . . . . . 619is.name : . 620is.nan : . . . . . . . 621is.null : . 622is.numeric : | 623is.numeric_version : . 624is.object : . 625is.ordered : . 626is.package_version : . 627is.pairlist : . 628is.primitive : . 629is.qr : . 630is.raw : . 631is.recursive : . 632is.symbol : . 633is.table : . 634is.vector : | 635 636> is.ALL(structure(1:7, names = paste("a",1:7,sep="")), true.only = TRUE) 637[1] "is.atomic" "is.integer" "is.numeric" "is.vector" 638 639 640 demo(recursion) 641 ---- ~~~~~~~~~ 642 643> # Copyright (C) 1997-2005 The R Core Team 644> 645> ## Adaptive integration: Venables and Ripley pp. 105-110 646> ## This is the basic integrator. 647> 648> area <- function(f, a, b, ..., fa = f(a, ...), fb = f(b, ...), limit 649+ = 10, eps = 1.e-5) 650+ { 651+ h <- b - a 652+ d <- (a + b)/2 653+ fd <- f(d, ...) 654+ a1 <- ((fa + fb) * h)/2 655+ a2 <- ((fa + 4 * fd + fb) * h)/6 656+ if(abs(a1 - a2) < eps) 657+ return(a2) 658+ if(limit == 0) { 659+ warning(paste("iteration limit reached near x = ", 660+ d)) 661+ return(a2) 662+ } 663+ area(f, a, d, ..., fa = fa, fb = fd, limit = limit - 1, 664+ eps = eps) + area(f, d, b, ..., fa = fd, fb = 665+ fb, limit = limit - 1, eps = eps) 666+ } 667 668> ## The function to be integrated 669> 670> fbeta <- function(x, alpha, beta) 671+ { 672+ x^(alpha - 1) * (1 - x)^(beta - 1) 673+ } 674 675> ## Compute the approximate integral, the exact integral and the error 676> 677> b0 <- area(fbeta, 0, 1, alpha=3.5, beta=1.5) 678 679> b1 <- exp(lgamma(3.5) + lgamma(1.5) - lgamma(5)) 680 681> c(b0, b1, b0-b1) 682[1] 1.227170e-01 1.227185e-01 -1.443996e-06 683 684> ## Modify the function so that it records where it was evaluated 685> 686> fbeta.tmp <- function (x, alpha, beta) 687+ { 688+ val <<- c(val, x) 689+ x^(alpha - 1) * (1 - x)^(beta - 1) 690+ } 691 692> ## Recompute and plot the evaluation points. 693> 694> val <- NULL 695 696> b0 <- area(fbeta.tmp, 0, 1, alpha=3.5, beta=1.5) 697 698> plot(val, fbeta(val, 3.5, 1.5), pch=0) 699 700> ## Better programming style -- renaming the function will have no effect. 701> ## The use of "Recall" as in V+R is VERY black magic. You can get the 702> ## same effect transparently by supplying a wrapper function. 703> ## This is the approved Abelson+Sussman method. 704> 705> area <- function(f, a, b, ..., limit=10, eps=1e-5) { 706+ area2 <- function(f, a, b, ..., fa = f(a, ...), fb = f(b, ...), 707+ limit = limit, eps = eps) { 708+ h <- b - a 709+ d <- (a + b)/2 710+ fd <- f(d, ...) 711+ a1 <- ((fa + fb) * h)/2 712+ a2 <- ((fa + 4 * fd + fb) * h)/6 713+ if(abs(a1 - a2) < eps) 714+ return(a2) 715+ if(limit == 0) { 716+ warning(paste("iteration limit reached near x =", d)) 717+ return(a2) 718+ } 719+ area2(f, a, d, ..., fa = fa, fb = fd, limit = limit - 1, 720+ eps = eps) + area2(f, d, b, ..., fa = fd, fb = 721+ fb, limit = limit - 1, eps = eps) 722+ } 723+ area2(f, a, b, ..., limit=limit, eps=eps) 724+ } 725 726 727 demo(scoping) 728 ---- ~~~~~~~ 729 730> ## Here is a little example which shows a fundamental difference between 731> ## R and S. It is a little example from Abelson and Sussman which models 732> ## the way in which bank accounts work. It shows how R functions can 733> ## encapsulate state information. 734> ## 735> ## When invoked, "open.account" defines and returns three functions 736> ## in a list. Because the variable "total" exists in the environment 737> ## where these functions are defined they have access to its value. 738> ## This is even true when "open.account" has returned. The only way 739> ## to access the value of "total" is through the accessor functions 740> ## withdraw, deposit and balance. Separate accounts maintain their 741> ## own balances. 742> ## 743> ## This is a very nifty way of creating "closures" and a little thought 744> ## will show you that there are many ways of using this in statistics. 745> 746> # Copyright (C) 1997-8 The R Core Team 747> 748> open.account <- function(total) { 749+ 750+ list( 751+ deposit = function(amount) { 752+ if(amount <= 0) 753+ stop("Deposits must be positive!\n") 754+ total <<- total + amount 755+ cat(amount,"deposited. Your balance is", total, "\n\n") 756+ }, 757+ withdraw = function(amount) { 758+ if(amount > total) 759+ stop("You don't have that much money!\n") 760+ total <<- total - amount 761+ cat(amount,"withdrawn. Your balance is", total, "\n\n") 762+ }, 763+ balance = function() { 764+ cat("Your balance is", total, "\n\n") 765+ } 766+ ) 767+ } 768 769> ross <- open.account(100) 770 771> robert <- open.account(200) 772 773> ross$withdraw(30) 77430 withdrawn. Your balance is 70 775 776 777> ross$balance() 778Your balance is 70 779 780 781> robert$balance() 782Your balance is 200 783 784 785> ross$deposit(50) 78650 deposited. Your balance is 120 787 788 789> ross$balance() 790Your balance is 120 791 792 793> try(ross$withdraw(500)) # no way.. 794Error in ross$withdraw(500) : You don't have that much money! 795 796 797 798 demo(graphics) 799 ---- ~~~~~~~~ 800 801> # Copyright (C) 1997-2009 The R Core Team 802> 803> require(datasets) 804 805> require(grDevices); require(graphics) 806 807> ## Here is some code which illustrates some of the differences between 808> ## R and S graphics capabilities. Note that colors are generally specified 809> ## by a character string name (taken from the X11 rgb.txt file) and that line 810> ## textures are given similarly. The parameter "bg" sets the background 811> ## parameter for the plot and there is also an "fg" parameter which sets 812> ## the foreground color. 813> 814> 815> x <- stats::rnorm(50) 816 817> opar <- par(bg = "white") 818 819> plot(x, ann = FALSE, type = "n") 820 821> abline(h = 0, col = gray(.90)) 822 823> lines(x, col = "green4", lty = "dotted") 824 825> points(x, bg = "limegreen", pch = 21) 826 827> title(main = "Simple Use of Color In a Plot", 828+ xlab = "Just a Whisper of a Label", 829+ col.main = "blue", col.lab = gray(.8), 830+ cex.main = 1.2, cex.lab = 1.0, font.main = 4, font.lab = 3) 831 832> ## A little color wheel. This code just plots equally spaced hues in 833> ## a pie chart. If you have a cheap SVGA monitor (like me) you will 834> ## probably find that numerically equispaced does not mean visually 835> ## equispaced. On my display at home, these colors tend to cluster at 836> ## the RGB primaries. On the other hand on the SGI Indy at work the 837> ## effect is near perfect. 838> 839> par(bg = "gray") 840 841> pie(rep(1,24), col = rainbow(24), radius = 0.9) 842 843> title(main = "A Sample Color Wheel", cex.main = 1.4, font.main = 3) 844 845> title(xlab = "(Use this as a test of monitor linearity)", 846+ cex.lab = 0.8, font.lab = 3) 847 848> ## We have already confessed to having these. This is just showing off X11 849> ## color names (and the example (from the postscript manual) is pretty "cute". 850> 851> pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12) 852 853> names(pie.sales) <- c("Blueberry", "Cherry", 854+ "Apple", "Boston Cream", "Other", "Vanilla Cream") 855 856> pie(pie.sales, 857+ col = c("purple","violetred1","green3","cornsilk","cyan","white")) 858 859> title(main = "January Pie Sales", cex.main = 1.8, font.main = 1) 860 861> title(xlab = "(Don't try this at home kids)", cex.lab = 0.8, font.lab = 3) 862 863> ## Boxplots: I couldn't resist the capability for filling the "box". 864> ## The use of color seems like a useful addition, it focuses attention 865> ## on the central bulk of the data. 866> 867> par(bg="cornsilk") 868 869> n <- 10 870 871> g <- gl(n, 100, n*100) 872 873> x <- rnorm(n*100) + sqrt(as.numeric(g)) 874 875> boxplot(split(x,g), col="lavender", notch=TRUE) 876 877> title(main="Notched Boxplots", xlab="Group", font.main=4, font.lab=1) 878 879> ## An example showing how to fill between curves. 880> 881> par(bg="white") 882 883> n <- 100 884 885> x <- c(0,cumsum(rnorm(n))) 886 887> y <- c(0,cumsum(rnorm(n))) 888 889> xx <- c(0:n, n:0) 890 891> yy <- c(x, rev(y)) 892 893> plot(xx, yy, type="n", xlab="Time", ylab="Distance") 894 895> polygon(xx, yy, col="gray") 896 897> title("Distance Between Brownian Motions") 898 899> ## Colored plot margins, axis labels and titles. You do need to be 900> ## careful with these kinds of effects. It's easy to go completely 901> ## over the top and you can end up with your lunch all over the keyboard. 902> ## On the other hand, my market research clients love it. 903> 904> x <- c(0.00, 0.40, 0.86, 0.85, 0.69, 0.48, 0.54, 1.09, 1.11, 1.73, 2.05, 2.02) 905 906> par(bg="lightgray") 907 908> plot(x, type="n", axes=FALSE, ann=FALSE) 909 910> usr <- par("usr") 911 912> rect(usr[1], usr[3], usr[2], usr[4], col="cornsilk", border="black") 913 914> lines(x, col="blue") 915 916> points(x, pch=21, bg="lightcyan", cex=1.25) 917 918> axis(2, col.axis="blue", las=1) 919 920> axis(1, at=1:12, lab=month.abb, col.axis="blue") 921 922> box() 923 924> title(main= "The Level of Interest in R", font.main=4, col.main="red") 925 926> title(xlab= "1996", col.lab="red") 927 928> ## A filled histogram, showing how to change the font used for the 929> ## main title without changing the other annotation. 930> 931> par(bg="cornsilk") 932 933> x <- rnorm(1000) 934 935> hist(x, xlim=range(-4, 4, x), col="lavender", main="") 936 937> title(main="1000 Normal Random Variates", font.main=3) 938 939> ## A scatterplot matrix 940> ## The good old Iris data (yet again) 941> 942> pairs(iris[1:4], main="Edgar Anderson's Iris Data", font.main=4, pch=19) 943 944> pairs(iris[1:4], main="Edgar Anderson's Iris Data", pch=21, 945+ bg = c("red", "green3", "blue")[unclass(iris$Species)]) 946 947> ## Contour plotting 948> ## This produces a topographic map of one of Auckland's many volcanic "peaks". 949> 950> x <- 10*1:nrow(volcano) 951 952> y <- 10*1:ncol(volcano) 953 954> lev <- pretty(range(volcano), 10) 955 956> par(bg = "lightcyan") 957 958> pin <- par("pin") 959 960> xdelta <- diff(range(x)) 961 962> ydelta <- diff(range(y)) 963 964> xscale <- pin[1]/xdelta 965 966> yscale <- pin[2]/ydelta 967 968> scale <- min(xscale, yscale) 969 970> xadd <- 0.5*(pin[1]/scale - xdelta) 971 972> yadd <- 0.5*(pin[2]/scale - ydelta) 973 974> plot(numeric(0), numeric(0), 975+ xlim = range(x)+c(-1,1)*xadd, ylim = range(y)+c(-1,1)*yadd, 976+ type = "n", ann = FALSE) 977 978> usr <- par("usr") 979 980> rect(usr[1], usr[3], usr[2], usr[4], col="green3") 981 982> contour(x, y, volcano, levels = lev, col="yellow", lty="solid", add=TRUE) 983 984> box() 985 986> title("A Topographic Map of Maunga Whau", font= 4) 987 988> title(xlab = "Meters North", ylab = "Meters West", font= 3) 989 990> mtext("10 Meter Contour Spacing", side=3, line=0.35, outer=FALSE, 991+ at = mean(par("usr")[1:2]), cex=0.7, font=3) 992 993> ## Conditioning plots 994> 995> par(bg="cornsilk") 996 997> coplot(lat ~ long | depth, data = quakes, pch = 21, bg = "green3") 998 999> par(opar) 1000 1001 1002 demo(image) 1003 ---- ~~~~~ 1004 1005> # Copyright (C) 1997-2009 The R Core Team 1006> 1007> require(datasets) 1008 1009> require(grDevices); require(graphics) 1010 1011> x <- 10*(1:nrow(volcano)); x.at <- seq(100, 800, by=100) 1012 1013> y <- 10*(1:ncol(volcano)); y.at <- seq(100, 600, by=100) 1014 1015> # Using Terrain Colors 1016> 1017> image(x, y, volcano, col=terrain.colors(100),axes=FALSE) 1018 1019> contour(x, y, volcano, levels=seq(90, 200, by=5), add=TRUE, col="brown") 1020 1021> axis(1, at=x.at) 1022 1023> axis(2, at=y.at) 1024 1025> box() 1026 1027> title(main="Maunga Whau Volcano", sub = "col=terrain.colors(100)", font.main=4) 1028 1029> # Using Heat Colors 1030> 1031> image(x, y, volcano, col=heat.colors(100), axes=FALSE) 1032 1033> contour(x, y, volcano, levels=seq(90, 200, by=5), add=TRUE, col="brown") 1034 1035> axis(1, at=x.at) 1036 1037> axis(2, at=y.at) 1038 1039> box() 1040 1041> title(main="Maunga Whau Volcano", sub = "col=heat.colors(100)", font.main=4) 1042 1043> # Using Gray Scale 1044> 1045> image(x, y, volcano, col=gray(100:200/200), axes=FALSE) 1046 1047> contour(x, y, volcano, levels=seq(90, 200, by=5), add=TRUE, col="black") 1048 1049> axis(1, at=x.at) 1050 1051> axis(2, at=y.at) 1052 1053> box() 1054 1055> title(main="Maunga Whau Volcano \n col=gray(100:200/200)", font.main=4) 1056 1057> ## Filled Contours are even nicer sometimes : 1058> example(filled.contour) 1059 1060flld.c> require("grDevices") # for colours 1061 1062flld.c> filled.contour(volcano, asp = 1) # simple 1063 1064flld.c> x <- 10*1:nrow(volcano) 1065 1066flld.c> y <- 10*1:ncol(volcano) 1067 1068flld.c> filled.contour(x, y, volcano, 1069flld.c+ color.palette = function(n) hcl.colors(n, "terrain"), 1070flld.c+ plot.title = title(main = "The Topography of Maunga Whau", 1071flld.c+ xlab = "Meters North", ylab = "Meters West"), 1072flld.c+ plot.axes = { axis(1, seq(100, 800, by = 100)) 1073flld.c+ axis(2, seq(100, 600, by = 100)) }, 1074flld.c+ key.title = title(main = "Height\n(meters)"), 1075flld.c+ key.axes = axis(4, seq(90, 190, by = 10))) # maybe also asp = 1 1076 1077flld.c> mtext(paste("filled.contour(.) from", R.version.string), 1078flld.c+ side = 1, line = 4, adj = 1, cex = .66) 1079 1080flld.c> # Annotating a filled contour plot 1081flld.c> a <- expand.grid(1:20, 1:20) 1082 1083flld.c> b <- matrix(a[,1] + a[,2], 20) 1084 1085flld.c> filled.contour(x = 1:20, y = 1:20, z = b, 1086flld.c+ plot.axes = { axis(1); axis(2); points(10, 10) }) 1087 1088flld.c> ## Persian Rug Art: 1089flld.c> x <- y <- seq(-4*pi, 4*pi, length.out = 27) 1090 1091flld.c> r <- sqrt(outer(x^2, y^2, "+")) 1092 1093flld.c> filled.contour(cos(r^2)*exp(-r/(2*pi)), axes = FALSE) 1094 1095flld.c> ## rather, the key *should* be labeled: 1096flld.c> filled.contour(cos(r^2)*exp(-r/(2*pi)), frame.plot = FALSE, 1097flld.c+ plot.axes = {}) 1098 1099 1100 demo(persp) 1101 ---- ~~~~~ 1102 1103> ### Demos for persp() plots -- things not in example(persp) 1104> ### ------------------------- 1105> 1106> require(datasets) 1107 1108> require(grDevices); require(graphics) 1109 1110> ## (1) The Obligatory Mathematical surface. 1111> ## Rotated sinc function. 1112> 1113> x <- seq(-10, 10, length.out = 50) 1114 1115> y <- x 1116 1117> rotsinc <- function(x,y) 1118+ { 1119+ sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y } 1120+ 10 * sinc( sqrt(x^2+y^2) ) 1121+ } 1122 1123> sinc.exp <- expression(z == Sinc(sqrt(x^2 + y^2))) 1124 1125> z <- outer(x, y, rotsinc) 1126 1127> oldpar <- par(bg = "white") 1128 1129> persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue") 1130 1131> title(sub=".")## work around persp+plotmath bug 1132 1133> title(main = sinc.exp) 1134 1135> persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue", 1136+ ltheta = 120, shade = 0.75, ticktype = "detailed", 1137+ xlab = "X", ylab = "Y", zlab = "Z") 1138 1139> title(sub=".")## work around persp+plotmath bug 1140 1141> title(main = sinc.exp) 1142 1143> ## (2) Visualizing a simple DEM model 1144> 1145> z <- 2 * volcano # Exaggerate the relief 1146 1147> x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N) 1148 1149> y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W) 1150 1151> persp(x, y, z, theta = 120, phi = 15, scale = FALSE, axes = FALSE) 1152 1153> ## (3) Now something more complex 1154> ## We border the surface, to make it more "slice like" 1155> ## and color the top and sides of the surface differently. 1156> 1157> z0 <- min(z) - 20 1158 1159> z <- rbind(z0, cbind(z0, z, z0), z0) 1160 1161> x <- c(min(x) - 1e-10, x, max(x) + 1e-10) 1162 1163> y <- c(min(y) - 1e-10, y, max(y) + 1e-10) 1164 1165> fill <- matrix("green3", nrow = nrow(z)-1, ncol = ncol(z)-1) 1166 1167> fill[ , i2 <- c(1,ncol(fill))] <- "gray" 1168 1169> fill[i1 <- c(1,nrow(fill)) , ] <- "gray" 1170 1171> par(bg = "lightblue") 1172 1173> persp(x, y, z, theta = 120, phi = 15, col = fill, scale = FALSE, axes = FALSE) 1174 1175> title(main = "Maunga Whau\nOne of 50 Volcanoes in the Auckland Region.", 1176+ font.main = 4) 1177 1178> par(bg = "slategray") 1179 1180> persp(x, y, z, theta = 135, phi = 30, col = fill, scale = FALSE, 1181+ ltheta = -120, lphi = 15, shade = 0.65, axes = FALSE) 1182 1183> ## Don't draw the grid lines : border = NA 1184> persp(x, y, z, theta = 135, phi = 30, col = "green3", scale = FALSE, 1185+ ltheta = -120, shade = 0.75, border = NA, box = FALSE) 1186 1187> ## `color gradient in the soil' : 1188> fcol <- fill ; fcol[] <- terrain.colors(nrow(fcol)) 1189 1190> persp(x, y, z, theta = 135, phi = 30, col = fcol, scale = FALSE, 1191+ ltheta = -120, shade = 0.3, border = NA, box = FALSE) 1192 1193> ## `image like' colors on top : 1194> fcol <- fill 1195 1196> zi <- volcano[ -1,-1] + volcano[ -1,-61] + 1197+ volcano[-87,-1] + volcano[-87,-61] ## / 4 1198 1199> fcol[-i1,-i2] <- 1200+ terrain.colors(20)[cut(zi, 1201+ stats::quantile(zi, seq(0,1, length.out = 21)), 1202+ include.lowest = TRUE)] 1203 1204> persp(x, y, 2*z, theta = 110, phi = 40, col = fcol, scale = FALSE, 1205+ ltheta = -120, shade = 0.4, border = NA, box = FALSE) 1206 1207> ## reset par(): 1208> par(oldpar) 1209 1210 1211 demo(glm.vr) 1212 ---- ~~~~~~ 1213 1214> # Copyright (C) 1997-2009 The R Core Team 1215> 1216> #### -*- R -*- 1217> require(stats) 1218 1219> Fr <- c(68,42,42,30, 37,52,24,43, 1220+ 66,50,33,23, 47,55,23,47, 1221+ 63,53,29,27, 57,49,19,29) 1222 1223> Temp <- gl(2, 2, 24, labels = c("Low", "High")) 1224 1225> Soft <- gl(3, 8, 24, labels = c("Hard","Medium","Soft")) 1226 1227> M.user <- gl(2, 4, 24, labels = c("N", "Y")) 1228 1229> Brand <- gl(2, 1, 24, labels = c("X", "M")) 1230 1231> detg <- data.frame(Fr,Temp, Soft,M.user, Brand) 1232 1233> detg.m0 <- glm(Fr ~ M.user*Temp*Soft + Brand, family = poisson, data = detg) 1234 1235> summary(detg.m0) 1236 1237Call: 1238glm(formula = Fr ~ M.user * Temp * Soft + Brand, family = poisson, 1239 data = detg) 1240 1241Deviance Residuals: 1242 Min 1Q Median 3Q Max 1243-2.20876 -0.99190 -0.00126 0.93542 1.97601 1244 1245Coefficients: 1246 Estimate Std. Error z value Pr(>|z|) 1247(Intercept) 4.01524 0.10034 40.018 < 2e-16 *** 1248M.userY -0.21184 0.14257 -1.486 0.13731 1249TempHigh -0.42381 0.15159 -2.796 0.00518 ** 1250SoftMedium 0.05311 0.13308 0.399 0.68984 1251SoftSoft 0.05311 0.13308 0.399 0.68984 1252BrandM -0.01587 0.06300 -0.252 0.80106 1253M.userY:TempHigh 0.13987 0.22168 0.631 0.52806 1254M.userY:SoftMedium 0.08323 0.19685 0.423 0.67245 1255M.userY:SoftSoft 0.12169 0.19591 0.621 0.53449 1256TempHigh:SoftMedium -0.30442 0.22239 -1.369 0.17104 1257TempHigh:SoftSoft -0.30442 0.22239 -1.369 0.17104 1258M.userY:TempHigh:SoftMedium 0.21189 0.31577 0.671 0.50220 1259M.userY:TempHigh:SoftSoft -0.20387 0.32540 -0.627 0.53098 1260--- 1261Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 1262 1263(Dispersion parameter for poisson family taken to be 1) 1264 1265 Null deviance: 118.627 on 23 degrees of freedom 1266Residual deviance: 32.826 on 11 degrees of freedom 1267AIC: 191.24 1268 1269Number of Fisher Scoring iterations: 4 1270 1271 1272> detg.mod <- glm(terms(Fr ~ M.user*Temp*Soft + Brand*M.user*Temp, 1273+ keep.order = TRUE), 1274+ family = poisson, data = detg) 1275 1276> summary(detg.mod) 1277 1278Call: 1279glm(formula = terms(Fr ~ M.user * Temp * Soft + Brand * M.user * 1280 Temp, keep.order = TRUE), family = poisson, data = detg) 1281 1282Deviance Residuals: 1283 Min 1Q Median 3Q Max 1284-0.91365 -0.35585 0.00253 0.33027 0.92146 1285 1286Coefficients: 1287 Estimate Std. Error z value Pr(>|z|) 1288(Intercept) 4.14887 0.10603 39.128 < 2e-16 *** 1289M.userY -0.40521 0.16188 -2.503 0.01231 * 1290TempHigh -0.44275 0.17121 -2.586 0.00971 ** 1291M.userY:TempHigh -0.12692 0.26257 -0.483 0.62883 1292SoftMedium 0.05311 0.13308 0.399 0.68984 1293SoftSoft 0.05311 0.13308 0.399 0.68984 1294M.userY:SoftMedium 0.08323 0.19685 0.423 0.67245 1295M.userY:SoftSoft 0.12169 0.19591 0.621 0.53449 1296TempHigh:SoftMedium -0.30442 0.22239 -1.369 0.17104 1297TempHigh:SoftSoft -0.30442 0.22239 -1.369 0.17104 1298M.userY:TempHigh:SoftMedium 0.21189 0.31577 0.671 0.50220 1299M.userY:TempHigh:SoftSoft -0.20387 0.32540 -0.627 0.53098 1300BrandM -0.30647 0.10942 -2.801 0.00510 ** 1301M.userY:BrandM 0.40757 0.15961 2.554 0.01066 * 1302TempHigh:BrandM 0.04411 0.18463 0.239 0.81119 1303M.userY:TempHigh:BrandM 0.44427 0.26673 1.666 0.09579 . 1304--- 1305Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 1306 1307(Dispersion parameter for poisson family taken to be 1) 1308 1309 Null deviance: 118.627 on 23 degrees of freedom 1310Residual deviance: 5.656 on 8 degrees of freedom 1311AIC: 170.07 1312 1313Number of Fisher Scoring iterations: 4 1314 1315 1316> summary(detg.mod, correlation = TRUE, symbolic.cor = TRUE) 1317 1318Call: 1319glm(formula = terms(Fr ~ M.user * Temp * Soft + Brand * M.user * 1320 Temp, keep.order = TRUE), family = poisson, data = detg) 1321 1322Deviance Residuals: 1323 Min 1Q Median 3Q Max 1324-0.91365 -0.35585 0.00253 0.33027 0.92146 1325 1326Coefficients: 1327 Estimate Std. Error z value Pr(>|z|) 1328(Intercept) 4.14887 0.10603 39.128 < 2e-16 *** 1329M.userY -0.40521 0.16188 -2.503 0.01231 * 1330TempHigh -0.44275 0.17121 -2.586 0.00971 ** 1331M.userY:TempHigh -0.12692 0.26257 -0.483 0.62883 1332SoftMedium 0.05311 0.13308 0.399 0.68984 1333SoftSoft 0.05311 0.13308 0.399 0.68984 1334M.userY:SoftMedium 0.08323 0.19685 0.423 0.67245 1335M.userY:SoftSoft 0.12169 0.19591 0.621 0.53449 1336TempHigh:SoftMedium -0.30442 0.22239 -1.369 0.17104 1337TempHigh:SoftSoft -0.30442 0.22239 -1.369 0.17104 1338M.userY:TempHigh:SoftMedium 0.21189 0.31577 0.671 0.50220 1339M.userY:TempHigh:SoftSoft -0.20387 0.32540 -0.627 0.53098 1340BrandM -0.30647 0.10942 -2.801 0.00510 ** 1341M.userY:BrandM 0.40757 0.15961 2.554 0.01066 * 1342TempHigh:BrandM 0.04411 0.18463 0.239 0.81119 1343M.userY:TempHigh:BrandM 0.44427 0.26673 1.666 0.09579 . 1344--- 1345Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 1346 1347(Dispersion parameter for poisson family taken to be 1) 1348 1349 Null deviance: 118.627 on 23 degrees of freedom 1350Residual deviance: 5.656 on 8 degrees of freedom 1351AIC: 170.07 1352 1353Number of Fisher Scoring iterations: 4 1354 1355Correlation of Coefficients: 1356 1357(Intercept) 1 1358M.userY , 1 1359TempHigh , . 1 1360M.userY:TempHigh . , , 1 1361SoftMedium , . . 1 1362SoftSoft , . . . 1 1363M.userY:SoftMedium . , . , . 1 1364M.userY:SoftSoft . , . . , . 1 1365TempHigh:SoftMedium . , . . . . 1 1366TempHigh:SoftSoft . , . . . . . 1 1367M.userY:TempHigh:SoftMedium . . . . , . , . 1 1368M.userY:TempHigh:SoftSoft . . . . . , . , . 1 1369BrandM . 1 1370M.userY:BrandM . , 1 1371TempHigh:BrandM . . . . 1 1372M.userY:TempHigh:BrandM . . . . , 1 1373attr(,"legend") 1374[1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1 1375 1376 1377> anova(detg.m0, detg.mod) 1378Analysis of Deviance Table 1379 1380Model 1: Fr ~ M.user * Temp * Soft + Brand 1381Model 2: Fr ~ M.user * Temp * Soft + Brand * M.user * Temp 1382 Resid. Df Resid. Dev Df Deviance 13831 11 32.826 13842 8 5.656 3 27.17 1385 1386 1387 demo(smooth) 1388 ---- ~~~~~~ 1389 1390> ### This used to be in example(smooth) before we had package-specific demos 1391> # Copyright (C) 1997-2009 The R Core Team 1392> 1393> require(stats); require(graphics); require(datasets) 1394 1395> op <- par(mfrow = c(1,1)) 1396 1397> ## The help(smooth) examples: 1398> example(smooth, package="stats") 1399 1400smooth> require(graphics) 1401 1402smooth> ## see also demo(smooth) ! 1403smooth> 1404smooth> x1 <- c(4, 1, 3, 6, 6, 4, 1, 6, 2, 4, 2) # very artificial 1405 1406smooth> (x3R <- smooth(x1, "3R")) # 2 iterations of "3" 14073R Tukey smoother resulting from smooth(x = x1, kind = "3R") 1408 used 2 iterations 1409 [1] 3 3 3 6 6 4 4 4 2 2 2 1410 1411smooth> smooth(x3R, kind = "S") 1412S Tukey smoother resulting from smooth(x = x3R, kind = "S") 1413 changed 1414 [1] 3 3 3 3 4 4 4 4 2 2 2 1415 1416smooth> sm.3RS <- function(x, ...) 1417smooth+ smooth(smooth(x, "3R", ...), "S", ...) 1418 1419smooth> y <- c(1, 1, 19:1) 1420 1421smooth> plot(y, main = "misbehaviour of \"3RSR\"", col.main = 3) 1422 1423smooth> lines(sm.3RS(y)) 1424 1425smooth> lines(smooth(y)) 1426 1427smooth> lines(smooth(y, "3RSR"), col = 3, lwd = 2) # the horror 1428 1429smooth> x <- c(8:10, 10, 0, 0, 9, 9) 1430 1431smooth> plot(x, main = "breakdown of 3R and S and hence 3RSS") 1432 1433smooth> matlines(cbind(smooth(x, "3R"), smooth(x, "S"), smooth(x, "3RSS"), smooth(x))) 1434 1435smooth> presidents[is.na(presidents)] <- 0 # silly 1436 1437smooth> summary(sm3 <- smooth(presidents, "3R")) 14383R Tukey smoother resulting from 1439 smooth(x = presidents, kind = "3R") ; n = 120 1440 used 4 iterations 1441 Min. 1st Qu. Median Mean 3rd Qu. Max. 1442 0.0 44.0 57.0 54.2 71.0 82.0 1443 1444smooth> summary(sm2 <- smooth(presidents,"3RSS")) 14453RSS Tukey smoother resulting from 1446 smooth(x = presidents, kind = "3RSS") ; n = 120 1447 used 5 iterations 1448 Min. 1st Qu. Median Mean 3rd Qu. Max. 1449 0.00 44.00 57.00 55.45 69.00 82.00 1450 1451smooth> summary(sm <- smooth(presidents)) 14523RS3R Tukey smoother resulting from 1453 smooth(x = presidents) ; n = 120 1454 used 7 iterations 1455 Min. 1st Qu. Median Mean 3rd Qu. Max. 1456 24.00 44.00 57.00 55.88 69.00 82.00 1457 1458smooth> all.equal(c(sm2), c(smooth(smooth(sm3, "S"), "S"))) # 3RSS === 3R S S 1459[1] TRUE 1460 1461smooth> all.equal(c(sm), c(smooth(smooth(sm3, "S"), "3R"))) # 3RS3R === 3R S 3R 1462[1] TRUE 1463 1464smooth> plot(presidents, main = "smooth(presidents0, *) : 3R and default 3RS3R") 1465 1466smooth> lines(sm3, col = 3, lwd = 1.5) 1467 1468smooth> lines(sm, col = 2, lwd = 1.25) 1469 1470> ## Didactical investigation: 1471> 1472> showSmooth <- function(x, leg.x = 1, leg.y = max(x)) { 1473+ ss <- cbind(x, "3c" = smooth(x, "3", end="copy"), 1474+ "3" = smooth(x, "3"), 1475+ "3Rc" = smooth(x, "3R", end="copy"), 1476+ "3R" = smooth(x, "3R"), 1477+ sm = smooth(x)) 1478+ k <- ncol(ss) - 1 1479+ n <- length(x) 1480+ slwd <- c(1,1,4,1,3,2) 1481+ slty <- c(0, 2:(k+1)) 1482+ matplot(ss, main = "Tukey Smoothers", ylab = "y ; sm(y)", 1483+ type= c("p",rep("l",k)), pch= par("pch"), lwd= slwd, lty= slty) 1484+ legend(leg.x, leg.y, 1485+ c("Data", "3 (copy)", "3 (Tukey)", 1486+ "3R (copy)", "3R (Tukey)", "smooth()"), 1487+ pch= c(par("pch"),rep(-1,k)), col=1:(k+1), lwd= slwd, lty= slty) 1488+ ss 1489+ } 1490 1491> ## 4 simple didactical examples, showing different steps in smooth(): 1492> 1493> for(x in list(c(4, 6, 2, 2, 6, 3, 6, 6, 5, 2), 1494+ c(3, 2, 1, 4, 5, 1, 3, 2, 4, 5, 2), 1495+ c(2, 4, 2, 6, 1, 1, 2, 6, 3, 1, 6), 1496+ x1)) 1497+ print(t(showSmooth(x))) 1498 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] 1499x 4 6 2 2 6 3 6 6 5 2 15003c 4 4 2 2 3 6 6 6 5 2 15013 4 4 2 2 3 6 6 6 5 3 15023Rc 4 4 2 2 3 6 6 6 5 2 15033R 4 4 2 2 3 6 6 6 5 3 1504sm 4 4 4 3 3 6 6 6 5 3 1505 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] 1506x 3 2 1 4 5 1 3 2 4 5 2 15073c 3 2 2 4 4 3 2 3 4 4 2 15083 2 2 2 4 4 3 2 3 4 4 4 15093Rc 3 2 2 4 4 3 3 3 4 4 2 15103R 2 2 2 4 4 3 3 3 4 4 4 1511sm 2 2 2 2 3 3 3 3 4 4 4 1512 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] 1513x 2 4 2 6 1 1 2 6 3 1 6 15143c 2 2 4 2 1 1 2 3 3 3 6 15153 2 2 4 2 1 1 2 3 3 3 3 15163Rc 2 2 2 2 1 1 2 3 3 3 6 15173R 2 2 2 2 1 1 2 3 3 3 3 1518sm 2 2 2 2 2 2 2 3 3 3 3 1519 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] 1520x 4 1 3 6 6 4 1 6 2 4 2 15213c 4 3 3 6 6 4 4 2 4 2 2 15223 3 3 3 6 6 4 4 2 4 2 2 15233Rc 4 3 3 6 6 4 4 4 2 2 2 15243R 3 3 3 6 6 4 4 4 2 2 2 1525sm 3 3 3 3 4 4 4 4 2 2 2 1526 1527> par(op) 1528> 1529> cat("Time elapsed: ", proc.time() - .ptime, "\n") 1530Time elapsed: 1.761 0.058 1.828 0 0 1531> 1532