R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) 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(-28205.97837 + ,-53575.97509 + ,-10734.10365 + ,-28203.33237 + ,-12632.63316 + ,-49816.6061 + ,128260.4167 + ,-24410.76654 + ,-13726.68326 + ,15932.90644 + ,-3370.04581 + ,-48240.11541 + ,-18597.64035 + ,170.9164267 + ,-593.6191667 + ,28581.11289 + ,-10862.81869 + ,56745.20595 + ,-43326.26972 + ,-28314.27258 + ,-10527.5927 + ,112650.5268 + ,3835.427347 + ,-67833.84014 + ,-30399.44013 + ,2129.846585 + ,27214.00211 + ,-2323.270518 + ,-378.5910857 + ,20926.1108 + ,20499.3727 + ,-40978.72116 + ,17743.14774 + ,-22915.10356 + ,37187.6772 + ,45371.96532 + ,85466.19006 + ,4922.152679 + ,21691.73737 + ,44977.99426 + ,66997.89787 + ,-23361.61903 + ,3531.592198 + ,-23771.30324 + ,2125.131416 + ,148668.6976 + ,-18965.08258 + ,-21764.15738 + ,29919.60747 + ,-36003.53677 + ,-15211.26432 + ,207.4763478 + ,46102.30746 + ,-15373.87645 + ,21286.84469 + ,29228.8055 + ,-12389.70578 + ,-40668.40972 + ,-4430.810792 + ,20596.02624 + ,49755.73325 + ,19229.70659 + ,10811.52545 + ,23291.1477 + ,-5366.322956 + ,-3.80539097 + ,-32629.69244 + ,-12502.94656 + ,-31685.13626 + ,-35310.30862 + ,-19361.64567 + ,-17022.24096 + ,-15481.29004 + ,-2988.376826 + ,-55745.84757 + ,27951.11689 + ,38171.43778 + ,-30073.95166 + ,-47615.53195 + ,-11612.6142 + ,3218.025406 + ,101842.9026 + ,26171.54851 + ,-43604.84933 + ,448.3626342 + ,11495.27867 + ,-59640.42044 + ,81817.93506 + ,43729.61588 + ,-74050.8687 + ,14335.2942 + ,-46944.19132 + ,-9615.056132 + ,21970.60282 + ,32397.97398 + ,-5936.788519 + ,-42154.70614 + ,-18050.88949 + ,-5920.334193 + ,7318.700408 + ,44.02752798 + ,-28349.88514 + ,22066.43591 + ,-61260.59918 + ,27731.54205 + ,-53561.10647 + ,-1597.746573 + ,15142.77311 + ,-47923.03942 + ,-21820.33291 + ,17663.79566 + ,-17911.80513 + ,4690.44137 + ,-16477.17596 + ,-131.6358949 + ,-39408.41099 + ,31845.05861 + ,9834.410265 + ,53462.98908 + ,-24998.75087 + ,33892.1052 + ,-18094.48821 + ,-24353.38737 + ,-30618.84295 + ,-4297.236256 + ,-25909.32575 + ,-10037.17545 + ,-59422.20557 + ,51539.52667 + ,-792.4269494 + ,-1244.156816 + ,7468.84447 + ,-14060.2665 + ,32263.31756 + ,53907.47149 + ,-30118.14216 + ,54202.49053 + ,-28578.95619 + ,2426.764874 + ,-31463.24786 + ,7944.278771 + ,101122.8747 + ,-53265.06565 + ,86069.11802 + ,28256.20668 + ,-80876.21826 + ,-14505.98303 + ,-39679.02793 + ,-8865.163215 + ,2040.919633 + ,-8167.690543 + ,-7810.690543 + ,-8332.298618 + ,-8265.690543 + ,4012.527536 + ,7221.684806 + ,-8265.690543 + ,-8062.690543 + ,-5547.8152 + ,14800.55903 + ,3538.264049 + ,-5340.524358 + ,-7296.690543 + ,11820.15553) > 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] 3.171344e-07 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 2931.519 > (armose <- arm / armse) [1] 1.081809e-10 > (geo <- geomean(x)) [1] NaN Warning message: In log(x) : NaNs produced > (har <- harmean(x)) [1] -677.6357 > (qua <- quamean(x)) [1] 37427.13 > (win <- winmean(x)) [,1] [,2] [1,] -82.82275 2888.3886 [2,] -197.36984 2827.5856 [3,] -274.82807 2763.9879 [4,] -252.87317 2754.7273 [5,] -705.17628 2655.9238 [6,] -592.73371 2633.7391 [7,] -655.83541 2591.1716 [8,] -1378.03875 2458.5895 [9,] -1924.44034 2365.5454 [10,] -1869.21204 2314.8611 [11,] -1783.25967 2298.6458 [12,] -1792.58209 2290.9512 [13,] -1920.67608 2265.6198 [14,] -2015.64156 2236.7081 [15,] -2044.36874 2153.1293 [16,] -2088.44313 2140.1792 [17,] -2007.83903 2120.3936 [18,] -2015.78465 2086.9778 [19,] -2623.76920 1999.1526 [20,] -2623.08368 1970.2011 [21,] -3010.42549 1915.0646 [22,] -2754.10631 1839.4857 [23,] -2675.77003 1826.5516 [24,] -2344.69336 1777.3312 [25,] -2494.22003 1728.1304 [26,] -2568.55999 1711.9453 [27,] -2536.17443 1685.5279 [28,] -2554.18720 1675.4332 [29,] -2558.39429 1664.2937 [30,] -2590.47679 1658.9617 [31,] -2405.71386 1619.6156 [32,] -2564.42265 1592.5093 [33,] -3136.84955 1528.9520 [34,] -3368.30198 1500.0063 [35,] -3388.18946 1497.8230 [36,] -2945.84140 1441.2124 [37,] -2831.75456 1411.6976 [38,] -2779.09898 1389.8771 [39,] -2843.94965 1380.3779 [40,] -2725.55194 1364.4153 [41,] -2940.54742 1321.8069 [42,] -3206.89999 1272.1277 [43,] -2940.66219 1243.0010 [44,] -3389.97562 1194.9286 [45,] -2947.55473 1112.3831 [46,] -2932.31049 1092.7758 [47,] -2960.34502 1070.0677 [48,] -3549.21912 985.0151 [49,] -3633.25954 974.5658 [50,] -3799.31712 950.9589 [51,] -3826.54335 897.3705 [52,] -4253.02785 826.5878 [53,] -4084.83312 783.7224 [54,] -4098.90291 776.1020 > (tri <- trimean(x)) [,1] [,2] [1,] -418.4721 2776.4518 [2,] -762.5127 2653.0138 [3,] -1055.8146 2553.0280 [4,] -1329.4937 2469.4985 [5,] -1616.1264 2380.7664 [6,] -1812.6999 2310.3798 [7,] -2035.0048 2238.3743 [8,] -2253.3290 2168.2939 [9,] -2376.2293 2116.7547 [10,] -2433.4001 2076.2687 [11,] -2498.5599 2039.9135 [12,] -2574.7347 2002.5902 [13,] -2652.1943 1962.8980 [14,] -2720.0500 1922.8751 [15,] -2781.6294 1882.7600 [16,] -2842.6955 1849.0092 [17,] -2902.1654 1813.5261 [18,] -2969.5686 1776.8021 [19,] -3038.5372 1740.1788 [20,] -3067.4089 1709.6502 [21,] -3097.2734 1679.0899 [22,] -3102.9254 1651.1693 [23,] -3124.9618 1627.9644 [24,] -3152.5733 1603.5307 [25,] -3200.9988 1581.1949 [26,] -3242.3959 1561.1405 [27,] -3281.0354 1540.3337 [28,] -3322.9274 1519.6039 [29,] -3365.4050 1497.3621 [30,] -3409.2876 1473.4983 [31,] -3453.1715 1447.2177 [32,] -3508.5854 1421.5095 [33,] -3557.9613 1395.2892 [34,] -3579.7613 1372.2177 [35,] -3590.6121 1348.9958 [36,] -3600.9218 1322.8047 [37,] -3634.0802 1298.8013 [38,] -3674.4922 1274.4117 [39,] -3719.4263 1248.7447 [40,] -3763.2535 1220.2737 [41,] -3815.1386 1189.1665 [42,] -3858.8682 1158.2998 [43,] -3891.5064 1128.5514 [44,] -3939.2232 1097.2985 [45,] -3966.8880 1067.0313 [46,] -4018.4839 1041.7825 [47,] -4018.4839 1014.4236 [48,] -4130.9407 984.9926 [49,] -4161.0551 962.0782 [50,] -4188.6566 936.0377 [51,] -4209.2539 908.1782 [52,] -4229.7652 882.6521 [53,] -4228.5003 862.5638 [54,] -4236.4388 844.4358 > (midr <- midrange(x)) [1] 33896.24 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] -4055.574 -3815.139 -4055.574 -3815.139 -3815.139 -4055.574 -3815.139 [8] -3763.254 > postscript(file="/var/wessaorg/rcomp/tmp/1ysg41321547877.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2k77d1321547877.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/3q4il1321547877.tab") > > try(system("convert tmp/1ysg41321547877.ps tmp/1ysg41321547877.png",intern=TRUE)) character(0) > try(system("convert tmp/2k77d1321547877.ps tmp/2k77d1321547877.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.005 0.169 1.169