R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree" Copyright (C) 2015 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-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(250785 + ,250140 + ,255755 + ,254671 + ,253919 + ,253741 + ,252729 + ,253810 + ,256653 + ,255231 + ,258405 + ,251061 + ,254811 + ,254895 + ,258325 + ,257608 + ,258759 + ,258621 + ,257852 + ,260560 + ,262358 + ,260812 + ,261165 + ,257164 + ,260720 + ,259581 + ,264743 + ,261845 + ,262262 + ,261631 + ,258953 + ,259966 + ,262850 + ,262204 + ,263418 + ,262752 + ,266433 + ,267722 + ,266003 + ,262971 + ,265521 + ,264676 + ,270223 + ,269508 + ,268457 + ,265814 + ,266680 + ,263018 + ,269285 + ,269829 + ,270911 + ,266844 + ,271244 + ,269907 + ,271296 + ,270157 + ,271322 + ,267179 + ,264101 + ,265518 + ,269419 + ,268714 + ,272482 + ,268351 + ,268175 + ,270674 + ,272764 + ,272599 + ,270333 + ,270846 + ,270491 + ,269160 + ,274027 + ,273784 + ,276663 + ,274525 + ,271344 + ,271115 + ,270798 + ,273911 + ,273985 + ,271917 + ,273338 + ,270601 + ,273547 + ,275363 + ,281229 + ,277793 + ,279913 + ,282500 + ,280041 + ,282166 + ,290304 + ,283519 + ,287816 + ,285226 + ,287595 + ,289741 + ,289148 + ,288301 + ,290155 + ,289648 + ,288225 + ,289351 + ,294735 + ,305333 + ,309030 + ,310215 + ,321935 + ,325734 + ,320846 + ,323023 + ,319753 + ,321753 + ,320757 + ,324479 + ,324641 + ,322767 + ,324181 + ,321389 + ,327897 + ,334287 + ,332653 + ,334819 + ,335264 + ,339622 + ,342440 + ,346585 + ,335378 + ,337010 + ,339130 + ,341193 + ,343507 + ,348915 + ,346431 + ,348322 + ,348288 + ,346597 + ,351076 + ,355215 + ,350562 + ,355266 + ,361565 + ,363462 + ,366183 + ,365423 + ,369208 + ,366713 + ,369354 + ,371970 + ,371824 + ,373187 + ,367270 + ,368140 + ,373742 + ,364815 + ,368558 + ,371503 + ,372611 + ,370197 + ,375441 + ,375888 + ,375132 + ,381142 + ,372024 + ,376070 + ,376864 + ,371401 + ,375687 + ,384304 + ,380738 + ,379908 + ,384007 + ,384499 + ,385106 + ,387935 + ,380435 + ,379281 + ,384153 + ,380599) > ylimmax = '' > ylimmin = '' > main = 'Robustness of Central Tendency' > #'GNU S' R Code compiled by R2WASP v. 1.2.291 () > #Author: root > #To cite this work: Wessa, P., (2012), Central Tendency (v1.0.4) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_centraltendency.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > # > 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] 304357.4 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 3388.88 > (armose <- arm / armse) [1] 89.81063 > (geo <- geomean(x)) [1] 301120.8 > (har <- harmean(x)) [1] 298054.1 > (qua <- quamean(x)) [1] 307716.1 > (win <- winmean(x)) [,1] [,2] [1,] 304345.3 3386.432 [2,] 304341.6 3385.266 [3,] 304366.2 3382.429 [4,] 304385.3 3380.086 [5,] 304383.2 3379.391 [6,] 304291.3 3366.713 [7,] 304304.8 3362.285 [8,] 304304.9 3360.989 [9,] 304300.9 3359.604 [10,] 304290.3 3354.378 [11,] 304284.0 3346.962 [12,] 304182.7 3322.211 [13,] 304162.3 3312.292 [14,] 304182.6 3307.849 [15,] 304186.2 3304.227 [16,] 304206.4 3298.305 [17,] 304184.8 3294.202 [18,] 304067.4 3275.906 [19,] 304023.4 3267.834 [20,] 303980.9 3258.619 [21,] 303985.7 3244.951 [22,] 304026.2 3240.592 [23,] 304083.4 3232.670 [24,] 304061.9 3226.063 [25,] 304060.5 3223.451 [26,] 303937.6 3199.426 [27,] 303881.1 3179.633 [28,] 303891.7 3174.555 [29,] 303844.8 3158.272 [30,] 303784.8 3149.597 [31,] 303651.5 3131.324 [32,] 303622.5 3114.955 [33,] 303543.3 3102.661 [34,] 303422.6 3084.831 [35,] 303313.5 3070.911 [36,] 303122.9 3034.889 [37,] 302873.4 2981.560 [38,] 301665.0 2829.152 [39,] 301668.4 2826.922 [40,] 300920.9 2718.278 [41,] 300804.5 2706.181 [42,] 300488.5 2661.936 [43,] 300392.0 2644.282 [44,] 300488.8 2635.828 [45,] 300127.8 2588.758 [46,] 300166.7 2585.429 [47,] 300213.9 2575.113 [48,] 299579.0 2487.174 [49,] 299411.9 2449.789 [50,] 299114.4 2412.434 [51,] 298699.3 2367.109 [52,] 298631.4 2348.106 [53,] 298138.5 2278.949 [54,] 297686.4 2229.837 [55,] 297692.5 2223.637 [56,] 297581.7 2208.616 [57,] 297514.9 2185.559 [58,] 297013.5 2134.475 [59,] 295536.6 1984.879 [60,] 294837.6 1918.021 > (tri <- trimean(x)) [,1] [,2] [1,] 304192.5 3380.770 [2,] 304036.2 3374.357 [3,] 303878.2 3367.767 [4,] 303708.0 3361.382 [5,] 303528.7 3354.813 [6,] 303345.6 3347.534 [7,] 303174.7 3341.875 [8,] 302997.5 3336.132 [9,] 302816.0 3329.719 [10,] 302630.4 3322.577 [11,] 302441.3 3315.136 [12,] 302248.0 3307.580 [13,] 302248.0 3301.859 [14,] 301868.0 3296.257 [15,] 301669.6 3290.124 [16,] 301465.5 3283.308 [17,] 301254.3 3275.942 [18,] 301038.8 3267.796 [19,] 300825.6 3260.109 [20,] 300609.2 3251.908 [21,] 300389.3 3243.213 [22,] 300162.6 3234.294 [23,] 299926.7 3224.288 [24,] 299680.3 3213.355 [25,] 299427.5 3201.319 [26,] 299427.5 3187.721 [27,] 298904.7 3174.185 [28,] 298637.2 3160.276 [29,] 298360.3 3144.703 [30,] 298076.7 3128.203 [31,] 297786.4 3110.053 [32,] 297492.8 3090.883 [33,] 297190.4 3070.305 [34,] 296881.0 3047.850 [35,] 296566.1 3023.772 [36,] 296244.8 2997.524 [37,] 295920.4 2970.739 [38,] 295595.2 2944.799 [39,] 295313.3 2928.781 [40,] 295020.0 2910.044 [41,] 294749.0 2897.667 [42,] 294472.1 2883.714 [43,] 294197.8 2870.813 [44,] 293915.9 2856.583 [45,] 293617.2 2839.801 [46,] 293321.2 2823.996 [47,] 293009.8 2804.867 [48,] 292681.3 2782.509 [49,] 292681.3 2764.611 [50,] 292042.3 2746.238 [51,] 291715.9 2727.462 [52,] 291715.9 2709.146 [53,] 291052.9 2688.067 [54,] 290718.7 2669.668 [55,] 290386.9 2652.094 [56,] 290035.3 2629.871 [57,] 289667.8 2603.316 [58,] 289280.6 2572.384 [59,] 288893.5 2540.762 [60,] 288555.8 2523.081 > (midr <- midrange(x)) [1] 319037.5 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 293318.4 293617.2 293318.4 293617.2 293617.2 293318.4 293617.2 293915.9 > postscript(file="/var/wessaorg/rcomp/tmp/1lppn1457477295.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/2c9rs1457477295.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/3q18y1457477295.tab") > > try(system("convert tmp/1lppn1457477295.ps tmp/1lppn1457477295.png",intern=TRUE)) character(0) > try(system("convert tmp/2c9rs1457477295.ps tmp/2c9rs1457477295.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.239 0.129 1.375