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(7612 + ,7381 + ,6978 + ,6819 + ,6688 + ,6454 + ,6679 + ,6921 + ,7807 + ,7898 + ,7832 + ,7384 + ,7620 + ,7281 + ,6929 + ,6587 + ,6071 + ,5928 + ,5964 + ,6374 + ,7160 + ,7213 + ,6890 + ,6525 + ,6739 + ,6580 + ,6391 + ,6254 + ,6114 + ,5978 + ,6315 + ,6427 + ,7132 + ,7292 + ,7708 + ,7525 + ,7450 + ,7526 + ,7263 + ,7070 + ,6893 + ,6781 + ,7188 + ,7015 + ,8273 + ,8470 + ,8230 + ,8137 + ,8122 + ,8367 + ,8141 + ,7750 + ,7504 + ,7330 + ,7608 + ,7647 + ,8942 + ,8865 + ,8320 + ,8207 + ,8105 + ,8290 + ,8162 + ,8051 + ,7699 + ,7440 + ,7656 + ,7549 + ,9086 + ,8942 + ,8764 + ,8500 + ,8239 + ,8443 + ,8349 + ,8288 + ,7970 + ,7496 + ,7745 + ,7543 + ,9036 + ,9075 + ,8859 + ,8605 + ,8419 + ,8495 + ,8284 + ,7582 + ,7691 + ,7046 + ,7442 + ,7596 + ,8597 + ,8436 + ,7881 + ,7477 + ,7508 + ,7361 + ,7299 + ,6914 + ,6768 + ,6746 + ,7052 + ,7139 + ,7714 + ,7750 + ,7622 + ,7424 + ,7444 + ,7208 + ,7128 + ,7022 + ,6688 + ,6199 + ,6400 + ,6474 + ,7182 + ,7330 + ,7410 + ,7442 + ,7753 + ,7762 + ,7814 + ,7838 + ,7298 + ,7155 + ,7076 + ,7450 + ,8216 + ,8246 + ,8335 + ,8171 + ,8485 + ,8435 + ,8369 + ,8210 + ,7888 + ,8061 + ,8139 + ,7837 + ,8943 + ,8523 + ,8104 + ,7969 + ,7921 + ,7930 + ,7706 + ,7552 + ,7379 + ,6946 + ,7128 + ,7393 + ,8092 + ,8004 + ,7903 + ,7710 + ,7867 + ,7860 + ,7723 + ,7477 + ,7126 + ,7161 + ,7162 + ,7406 + ,7944 + ,8084 + ,8088 + ,7972 + ,8184 + ,7914 + ,7845 + ,7610 + ,7278 + ,6883 + ,7123 + ,7182 + ,7912 + ,7893 + ,7671 + ,7403 + ,7663 + ,7589 + ,7450 + ,7069 + ,6670 + ,6285 + ,6506 + ,6539 + ,7291 + ,7391 + ,7126 + ,6752 + ,6835 + ,6664 + ,6562 + ,6174 + ,5741 + ,5398 + ,5203 + ,5673 + ,6379 + ,6418 + ,6272 + ,6059) > 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] 7464.848 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 52.82014 > (armose <- arm / armse) [1] 141.3258 > (geo <- geomean(x)) [1] 7425.659 > (har <- harmean(x)) [1] 7385.074 > (qua <- quamean(x)) [1] 7502.687 > (win <- winmean(x)) [,1] [,2] [1,] 7465.750 52.61865 [2,] 7468.064 52.07144 [3,] 7467.696 51.70671 [4,] 7471.343 51.12959 [5,] 7472.225 50.99971 [6,] 7470.373 50.62537 [7,] 7472.946 50.20275 [8,] 7469.691 49.64479 [9,] 7464.574 48.52608 [10,] 7467.123 48.08421 [11,] 7464.480 47.45614 [12,] 7466.363 46.88978 [13,] 7467.191 46.70987 [14,] 7467.397 46.52362 [15,] 7468.500 46.13127 [16,] 7471.010 45.34541 [17,] 7470.843 45.23436 [18,] 7471.814 45.09963 [19,] 7471.162 44.84487 [20,] 7468.025 44.13590 [21,] 7468.745 44.00689 [22,] 7469.716 43.47525 [23,] 7470.392 43.06028 [24,] 7472.392 42.46060 [25,] 7471.044 41.84237 [26,] 7472.574 41.61983 [27,] 7475.088 41.23448 [28,] 7476.049 40.82059 [29,] 7473.206 40.34730 [30,] 7483.500 39.05136 [31,] 7483.044 38.82723 [32,] 7482.260 38.47460 [33,] 7482.745 38.23413 [34,] 7482.245 38.18733 [35,] 7487.049 36.93451 [36,] 7485.990 36.59898 [37,] 7485.446 36.34058 [38,] 7484.515 35.68841 [39,] 7486.618 35.40885 [40,] 7493.676 34.64997 [41,] 7493.877 34.06779 [42,] 7500.260 32.82124 [43,] 7501.525 32.66571 [44,] 7499.583 32.37177 [45,] 7503.333 31.86730 [46,] 7504.010 31.64229 [47,] 7500.554 30.99971 [48,] 7502.201 30.42866 [49,] 7498.598 28.74901 [50,] 7499.824 27.27308 [51,] 7501.074 27.07754 [52,] 7506.936 26.52727 [53,] 7502.000 25.84134 [54,] 7502.794 25.14589 [55,] 7500.637 24.92060 [56,] 7500.363 24.62123 [57,] 7512.936 23.47955 [58,] 7511.230 23.19621 [59,] 7509.784 23.07613 [60,] 7508.902 22.90627 [61,] 7507.407 22.78300 [62,] 7506.495 22.50859 [63,] 7504.333 21.97858 [64,] 7507.157 21.39102 [65,] 7503.971 20.87525 [66,] 7502.029 20.66717 [67,] 7502.029 20.61418 [68,] 7507.029 19.94205 > (tri <- trimean(x)) [,1] [,2] [1,] 7468.020 51.52519 [2,] 7470.335 50.35641 [3,] 7471.505 49.41372 [4,] 7472.827 48.54807 [5,] 7473.216 47.79646 [6,] 7473.427 47.02691 [7,] 7473.974 46.28451 [8,] 7474.133 45.57066 [9,] 7474.742 44.90144 [10,] 7475.995 44.35527 [11,] 7476.989 43.83349 [12,] 7478.278 43.35633 [13,] 7479.416 42.91453 [14,] 7480.506 42.46227 [15,] 7481.603 41.99817 [16,] 7482.640 41.54193 [17,] 7483.512 41.13163 [18,] 7484.417 40.70172 [19,] 7485.277 40.25254 [20,] 7486.201 39.79283 [21,] 7487.346 39.36109 [22,] 7488.475 38.90724 [23,] 7489.576 38.46378 [24,] 7490.667 38.02030 [25,] 7491.675 37.59136 [26,] 7492.783 37.17765 [27,] 7493.840 36.74783 [28,] 7494.797 36.31331 [29,] 7495.733 35.87499 [30,] 7496.833 35.43576 [31,] 7496.833 35.06117 [32,] 7498.150 34.67021 [33,] 7498.884 34.27090 [34,] 7499.618 33.85354 [35,] 7500.396 33.40084 [36,] 7500.985 33.00388 [37,] 7501.638 32.59447 [38,] 7502.336 32.16509 [39,] 7503.095 31.74415 [40,] 7503.790 31.30319 [41,] 7504.213 30.87857 [42,] 7504.642 30.45619 [43,] 7504.822 30.08867 [44,] 7504.957 29.69351 [45,] 7505.175 29.27847 [46,] 7505.250 28.85866 [47,] 7505.300 28.40983 [48,] 7505.491 27.96369 [49,] 7505.623 27.51428 [50,] 7505.904 27.15512 [51,] 7506.147 26.87503 [52,] 7506.350 26.57415 [53,] 7506.327 26.28012 [54,] 7506.500 26.00458 [55,] 7506.649 25.75009 [56,] 7506.891 25.47724 [57,] 7507.156 25.19033 [58,] 7506.920 24.96416 [59,] 7506.744 24.72643 [60,] 7506.619 24.45991 [61,] 7506.524 24.16512 [62,] 7506.524 23.83292 [63,] 7506.487 23.47340 [64,] 7506.579 23.11031 [65,] 7506.554 22.74962 [66,] 7506.667 22.38312 [67,] 7506.871 21.97242 [68,] 7507.088 21.48941 > (midr <- midrange(x)) [1] 7144.5 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 7501.379 7506.147 7501.379 7506.147 7506.147 7501.379 7506.147 7505.904 > postscript(file="/var/wessaorg/rcomp/tmp/1n4cf1457197675.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/2g8sp1457197675.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/36z951457197675.tab") > > try(system("convert tmp/1n4cf1457197675.ps tmp/1n4cf1457197675.png",intern=TRUE)) character(0) > try(system("convert tmp/2g8sp1457197675.ps tmp/2g8sp1457197675.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.476 0.159 1.634