R version 2.9.0 (2009-04-17) Copyright (C) 2009 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. > x <- c(7.3,7.2,7.1,7,6.8,6.6,6.3,6.1,6,6.4,6.8,6.8,6.7,6.6,6.5,6.5,6.4,6.3,6.1,5.9,5.9,6.3,6.7,6.8,6.8,6.8,6.8,6.9,6.8,6.7,6.6,6.5,6.5,7,7.5,7.6,7.6,7.6,7.8,8,8,8,7.9,7.9,8,8.5,9.2,9.4,9.5,9.5,9.6,9.7,9.7,9.6,9.5,9.4,9.3,9.6,10.2,10.2,10.1,9.9,9.8,9.8,9.7,9.5,9.3,9.1,9,9.5,10,10.2,10.1,10,9.9,10,9.9,9.7,9.5,9.2,9,9.3,9.8,9.8,9.6,9.4,9.3,9.2,9.2,9,8.8,8.7,8.7,9.1,9.7,9.8,9.6,9.4,9.4,9.5,9.4,9.3,9.2,9,8.9,9.2,9.8,9.9,9.6,9.2,9.1,9.1,9,8.9,8.7,8.5,8.3,8.5,8.7,8.4,8.1,7.8,7.7,7.5,7.2,6.8,6.7,6.4,6.3,6.8,7.3,7.1,7,6.8,6.6,6.3,6.1,6.1,6.3,6.3,6,6.2,6.4,6.8,7.5,7.5,7.6,7.6,7.4,7.3,7.1,6.9,6.8,7.5,7.6,7.8,8,8.1,8.2,8.3,8.2,8,7.9,7.6,7.6,8.3,8.4,8.4,8.4,8.4,8.6,8.9,8.8,8.3,7.5,7.2,7.4,8.8,9.3,9.3,8.7,8.2,8.3,8.5,8.6,8.5,8.2,8.1,7.9,8.6,8.7,8.7,8.5,8.4,8.5,8.7,8.7,8.6,8.5,8.3,8,8.2,8.1,8.1,8,7.9,7.9,8,8,7.9,8,7.7,7.2,7.5,7.3,7,7,7,7.2,7.3,7.1,6.8,6.4,6.1,6.5,7.7,7.9,7.5,6.9,6.6,6.9,7.7,8,8,7.7,7.3,7.4,8.1,8.3,8.2) > ylimmax = '' > ylimmin = '' > main = 'Robustness of Central Tendency' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Prof. Dr. P. Wessa > #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > #Technical description: Write here your technical program description (don't use hard returns!) > geomean <- function(x) { + return(exp(mean(log(x)))) + } > harmean <- function(x) { + return(1/mean(1/x)) + } > quamean <- function(x) { + return(sqrt(mean(x*x))) + } > winmean <- function(x) { + x <-sort(x[!is.na(x)]) + n<-length(x) + denom <- 3 + nodenom <- n/denom + if (nodenom>40) denom <- n/40 + sqrtn = sqrt(n) + roundnodenom = floor(nodenom) + win <- array(NA,dim=c(roundnodenom,2)) + for (j in 1:roundnodenom) { + win[j,1] <- (j*x[j+1]+sum(x[(j+1):(n-j)])+j*x[n-j])/n + win[j,2] <- sd(c(rep(x[j+1],j),x[(j+1):(n-j)],rep(x[n-j],j)))/sqrtn + } + return(win) + } > trimean <- function(x) { + x <-sort(x[!is.na(x)]) + n<-length(x) + denom <- 3 + nodenom <- n/denom + if (nodenom>40) denom <- n/40 + sqrtn = sqrt(n) + roundnodenom = floor(nodenom) + tri <- array(NA,dim=c(roundnodenom,2)) + for (j in 1:roundnodenom) { + tri[j,1] <- mean(x,trim=j/n) + tri[j,2] <- sd(x[(j+1):(n-j)]) / sqrt(n-j*2) + } + return(tri) + } > midrange <- function(x) { + return((max(x)+min(x))/2) + } > q1 <- function(data,n,p,i,f) { + np <- n*p; + i <<- floor(np) + f <<- np - i + qvalue <- (1-f)*data[i] + f*data[i+1] + } > q2 <- function(data,n,p,i,f) { + np <- (n+1)*p + i <<- floor(np) + f <<- np - i + qvalue <- (1-f)*data[i] + f*data[i+1] + } > q3 <- function(data,n,p,i,f) { + np <- n*p + i <<- floor(np) + f <<- np - i + if (f==0) { + qvalue <- data[i] + } else { + qvalue <- data[i+1] + } + } > q4 <- function(data,n,p,i,f) { + np <- n*p + i <<- floor(np) + f <<- np - i + if (f==0) { + qvalue <- (data[i]+data[i+1])/2 + } else { + qvalue <- data[i+1] + } + } > q5 <- function(data,n,p,i,f) { + np <- (n-1)*p + i <<- floor(np) + f <<- np - i + if (f==0) { + qvalue <- data[i+1] + } else { + qvalue <- data[i+1] + f*(data[i+2]-data[i+1]) + } + } > q6 <- function(data,n,p,i,f) { + np <- n*p+0.5 + i <<- floor(np) + f <<- np - i + qvalue <- data[i] + } > q7 <- function(data,n,p,i,f) { + np <- (n+1)*p + i <<- floor(np) + f <<- np - i + if (f==0) { + qvalue <- data[i] + } else { + qvalue <- f*data[i] + (1-f)*data[i+1] + } + } > q8 <- function(data,n,p,i,f) { + np <- (n+1)*p + i <<- floor(np) + f <<- np - i + if (f==0) { + qvalue <- data[i] + } else { + if (f == 0.5) { + qvalue <- (data[i]+data[i+1])/2 + } else { + if (f < 0.5) { + qvalue <- data[i] + } else { + qvalue <- data[i+1] + } + } + } + } > midmean <- function(x,def) { + x <-sort(x[!is.na(x)]) + n<-length(x) + if (def==1) { + qvalue1 <- q1(x,n,0.25,i,f) + qvalue3 <- q1(x,n,0.75,i,f) + } + if (def==2) { + qvalue1 <- q2(x,n,0.25,i,f) + qvalue3 <- q2(x,n,0.75,i,f) + } + if (def==3) { + qvalue1 <- q3(x,n,0.25,i,f) + qvalue3 <- q3(x,n,0.75,i,f) + } + if (def==4) { + qvalue1 <- q4(x,n,0.25,i,f) + qvalue3 <- q4(x,n,0.75,i,f) + } + if (def==5) { + qvalue1 <- q5(x,n,0.25,i,f) + qvalue3 <- q5(x,n,0.75,i,f) + } + if (def==6) { + qvalue1 <- q6(x,n,0.25,i,f) + qvalue3 <- q6(x,n,0.75,i,f) + } + if (def==7) { + qvalue1 <- q7(x,n,0.25,i,f) + qvalue3 <- q7(x,n,0.75,i,f) + } + if (def==8) { + qvalue1 <- q8(x,n,0.25,i,f) + qvalue3 <- q8(x,n,0.75,i,f) + } + midm <- 0 + myn <- 0 + roundno4 <- round(n/4) + round3no4 <- round(3*n/4) + for (i in 1:n) { + if ((x[i]>=qvalue1) & (x[i]<=qvalue3)){ + midm = midm + x[i] + myn = myn + 1 + } + } + midm = midm / myn + return(midm) + } > (arm <- mean(x)) [1] 8.055 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 0.07422237 > (armose <- arm / armse) [1] 108.5252 > (geo <- geomean(x)) [1] 7.972207 > (har <- harmean(x)) [1] 7.888826 > (qua <- quamean(x)) [1] 8.136318 > (win <- winmean(x)) [,1] [,2] [1,] 8.055000 0.07422237 [2,] 8.055833 0.07412340 [3,] 8.054583 0.07397545 [4,] 8.056250 0.07378616 [5,] 8.054167 0.07355012 [6,] 8.054167 0.07355012 [7,] 8.054167 0.07355012 [8,] 8.050833 0.07318942 [9,] 8.054583 0.07278038 [10,] 8.058750 0.07234631 [11,] 8.058750 0.07234631 [12,] 8.053750 0.07182574 [13,] 8.053750 0.07182574 [14,] 8.053750 0.07182574 [15,] 8.053750 0.07182574 [16,] 8.053750 0.07182574 [17,] 8.060833 0.07111777 [18,] 8.053333 0.07036681 [19,] 8.053333 0.07036681 [20,] 8.053333 0.07036681 [21,] 8.053333 0.07036681 [22,] 8.062500 0.06948487 [23,] 8.052917 0.06855985 [24,] 8.052917 0.06855985 [25,] 8.052917 0.06855985 [26,] 8.052917 0.06855985 [27,] 8.064167 0.06751618 [28,] 8.064167 0.06751618 [29,] 8.052083 0.06638963 [30,] 8.052083 0.06638963 [31,] 8.052083 0.06638963 [32,] 8.065417 0.06519509 [33,] 8.065417 0.06519509 [34,] 8.065417 0.06519509 [35,] 8.065417 0.06519509 [36,] 8.065417 0.06254189 [37,] 8.065417 0.06254189 [38,] 8.065417 0.06254189 [39,] 8.065417 0.06254189 [40,] 8.065417 0.06254189 [41,] 8.065417 0.06254189 [42,] 8.047917 0.06100892 [43,] 8.047917 0.06100892 [44,] 8.047917 0.06100892 [45,] 8.047917 0.06100892 [46,] 8.047917 0.06100892 [47,] 8.047917 0.06100892 [48,] 8.047917 0.06100892 [49,] 8.027500 0.05928716 [50,] 8.048333 0.05751407 [51,] 8.048333 0.05751407 [52,] 8.048333 0.05751407 [53,] 8.048333 0.05751407 [54,] 8.070833 0.05566822 [55,] 8.070833 0.05566822 [56,] 8.047500 0.05372110 [57,] 8.047500 0.05372110 [58,] 8.047500 0.05372110 [59,] 8.047500 0.05372110 [60,] 8.047500 0.04967441 [61,] 8.047500 0.04967441 [62,] 8.047500 0.04967441 [63,] 8.047500 0.04967441 [64,] 8.074167 0.04758459 [65,] 8.047083 0.04541724 [66,] 8.047083 0.04541724 [67,] 8.047083 0.04541724 [68,] 8.018750 0.04323189 [69,] 8.047500 0.04099497 [70,] 8.047500 0.04099497 [71,] 8.017917 0.03876874 [72,] 8.017917 0.03876874 [73,] 8.017917 0.03876874 [74,] 8.017917 0.03876874 [75,] 8.049167 0.03639045 [76,] 8.049167 0.03639045 [77,] 8.049167 0.03639045 [78,] 8.081667 0.03401312 [79,] 8.081667 0.03401312 [80,] 8.048333 0.03152320 > (tri <- trimean(x)) [,1] [,2] [1,] 8.055042 0.07374437 [2,] 8.055085 0.07324136 [3,] 8.054701 0.07276497 [4,] 8.054741 0.07231727 [5,] 8.054348 0.07189832 [6,] 8.054386 0.07151027 [7,] 8.054425 0.07110044 [8,] 8.054464 0.07066757 [9,] 8.054955 0.07026531 [10,] 8.055000 0.06989652 [11,] 8.054587 0.06956090 [12,] 8.054167 0.06920495 [13,] 8.054206 0.06888470 [14,] 8.054245 0.06854444 [15,] 8.054286 0.06818300 [16,] 8.054327 0.06779907 [17,] 8.054369 0.06739127 [18,] 8.053922 0.06701953 [19,] 8.053960 0.06668671 [20,] 8.054000 0.06633189 [21,] 8.054040 0.06595369 [22,] 8.054082 0.06555059 [23,] 8.053608 0.06518669 [24,] 8.053646 0.06486525 [25,] 8.053684 0.06452112 [26,] 8.053723 0.06415280 [27,] 8.053763 0.06375866 [28,] 8.053261 0.06340736 [29,] 8.052747 0.06303055 [30,] 8.052778 0.06270061 [31,] 8.052809 0.06234575 [32,] 8.052841 0.06196422 [33,] 8.052299 0.06162972 [34,] 8.051744 0.06126901 [35,] 8.051176 0.06088018 [36,] 8.050595 0.06046109 [37,] 8.050000 0.06017347 [38,] 8.049390 0.05986104 [39,] 8.048765 0.05952189 [40,] 8.048125 0.05915397 [41,] 8.047468 0.05875498 [42,] 8.046795 0.05832243 [43,] 8.046753 0.05794905 [44,] 8.046711 0.05754291 [45,] 8.046667 0.05710124 [46,] 8.046622 0.05662097 [47,] 8.046575 0.05609867 [48,] 8.046528 0.05553051 [49,] 8.046479 0.05491216 [50,] 8.047143 0.05435261 [51,] 8.047101 0.05385972 [52,] 8.047059 0.05332048 [53,] 8.047015 0.05273031 [54,] 8.046970 0.05208402 [55,] 8.046154 0.05150280 [56,] 8.045312 0.05086432 [57,] 8.045238 0.05030204 [58,] 8.045161 0.04968237 [59,] 8.045082 0.04899895 [60,] 8.045000 0.04824439 [61,] 8.044915 0.04771465 [62,] 8.044828 0.04712650 [63,] 8.044737 0.04647306 [64,] 8.044643 0.04574634 [65,] 8.043636 0.04510256 [66,] 8.043519 0.04455994 [67,] 8.043396 0.04395168 [68,] 8.043269 0.04326923 [69,] 8.044118 0.04268629 [70,] 8.044000 0.04221924 [71,] 8.043878 0.04168920 [72,] 8.044792 0.04127544 [73,] 8.045745 0.04080142 [74,] 8.046739 0.04025869 [75,] 8.047778 0.03963717 [76,] 8.047727 0.03915838 [77,] 8.047674 0.03860519 [78,] 8.047619 0.03796593 [79,] 8.046341 0.03746704 [80,] 8.045000 0.03688418 > (midr <- midrange(x)) [1] 8.05 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 7.995238 8.045000 7.995238 8.045000 8.045000 7.995238 8.045000 8.029231 > postscript(file="/var/www/html/rcomp/tmp/1bt7g1262273221.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > lb <- win[,1] - 2*win[,2] > ub <- win[,1] + 2*win[,2] > if ((ylimmin == '') | (ylimmax == '')) plot(win[,1],type='b',main=main, xlab='j', pch=19, ylab='Winsorized Mean(j/n)', ylim=c(min(lb),max(ub))) else plot(win[,1],type='l',main=main, xlab='j', pch=19, ylab='Winsorized Mean(j/n)', ylim=c(ylimmin,ylimmax)) > lines(ub,lty=3) > lines(lb,lty=3) > grid() > dev.off() null device 1 > postscript(file="/var/www/html/rcomp/tmp/2tg5f1262273221.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > lb <- tri[,1] - 2*tri[,2] > ub <- tri[,1] + 2*tri[,2] > if ((ylimmin == '') | (ylimmax == '')) plot(tri[,1],type='b',main=main, xlab='j', pch=19, ylab='Trimmed Mean(j/n)', ylim=c(min(lb),max(ub))) else plot(tri[,1],type='l',main=main, xlab='j', pch=19, ylab='Trimmed Mean(j/n)', ylim=c(ylimmin,ylimmax)) > lines(ub,lty=3) > lines(lb,lty=3) > grid() > dev.off() null device 1 > > #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/rcomp/createtable") > > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Central Tendency - Ungrouped Data',4,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Measure',header=TRUE) > a<-table.element(a,'Value',header=TRUE) > a<-table.element(a,'S.E.',header=TRUE) > a<-table.element(a,'Value/S.E.',header=TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink('http://www.xycoon.com/arithmetic_mean.htm', 'Arithmetic Mean', 'click to view the definition of the Arithmetic Mean'),header=TRUE) > a<-table.element(a,arm) > a<-table.element(a,hyperlink('http://www.xycoon.com/arithmetic_mean_standard_error.htm', armse, 'click to view the definition of the Standard Error of the Arithmetic Mean')) > a<-table.element(a,armose) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink('http://www.xycoon.com/geometric_mean.htm', 'Geometric Mean', 'click to view the definition of the Geometric Mean'),header=TRUE) > a<-table.element(a,geo) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink('http://www.xycoon.com/harmonic_mean.htm', 'Harmonic Mean', 'click to view the definition of the Harmonic Mean'),header=TRUE) > a<-table.element(a,har) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink('http://www.xycoon.com/quadratic_mean.htm', 'Quadratic Mean', 'click to view the definition of the Quadratic Mean'),header=TRUE) > a<-table.element(a,qua) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > for (j in 1:length(win[,1])) { + a<-table.row.start(a) + mylabel <- paste('Winsorized Mean (',j) + mylabel <- paste(mylabel,'/') + mylabel <- paste(mylabel,length(win[,1])) + mylabel <- paste(mylabel,')') + a<-table.element(a,hyperlink('http://www.xycoon.com/winsorized_mean.htm', mylabel, 'click to view the definition of the Winsorized Mean'),header=TRUE) + a<-table.element(a,win[j,1]) + a<-table.element(a,win[j,2]) + a<-table.element(a,win[j,1]/win[j,2]) + a<-table.row.end(a) + } > for (j in 1:length(tri[,1])) { + a<-table.row.start(a) + mylabel <- paste('Trimmed Mean (',j) + mylabel <- paste(mylabel,'/') + mylabel <- paste(mylabel,length(tri[,1])) + mylabel <- paste(mylabel,')') + a<-table.element(a,hyperlink('http://www.xycoon.com/arithmetic_mean.htm', mylabel, 'click to view the definition of the Trimmed Mean'),header=TRUE) + a<-table.element(a,tri[j,1]) + a<-table.element(a,tri[j,2]) + a<-table.element(a,tri[j,1]/tri[j,2]) + a<-table.row.end(a) + } > a<-table.row.start(a) > a<-table.element(a,hyperlink('http://www.xycoon.com/median_1.htm', 'Median', 'click to view the definition of the Median'),header=TRUE) > a<-table.element(a,median(x)) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,hyperlink('http://www.xycoon.com/midrange.htm', 'Midrange', 'click to view the definition of the Midrange'),header=TRUE) > a<-table.element(a,midr) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_1.htm','Weighted Average at Xnp',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[1]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_2.htm','Weighted Average at X(n+1)p',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[2]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_3.htm','Empirical Distribution Function',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[3]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_4.htm','Empirical Distribution Function - Averaging',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[4]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_5.htm','Empirical Distribution Function - Interpolation',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[5]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_6.htm','Closest Observation',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[6]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_7.htm','True Basic - Statistics Graphics Toolkit',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[7]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > mymid <- hyperlink('http://www.xycoon.com/midmean.htm', 'Midmean', 'click to view the definition of the Midmean') > mylabel <- paste(mymid,hyperlink('http://www.xycoon.com/method_8.htm','MS Excel (old versions)',''),sep=' - ') > a<-table.element(a,mylabel,header=TRUE) > a<-table.element(a,midm[8]) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Number of observations',header=TRUE) > a<-table.element(a,length(x)) > a<-table.element(a,'') > a<-table.element(a,'') > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/rcomp/tmp/3u8ap1262273221.tab") > > try(system("convert tmp/1bt7g1262273221.ps tmp/1bt7g1262273221.png",intern=TRUE)) character(0) > try(system("convert tmp/2tg5f1262273221.ps tmp/2tg5f1262273221.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.265 0.361 1.712