R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 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(350840 + ,354950 + ,342750 + ,318610 + ,303890 + ,293480 + ,286820 + ,283650 + ,276190 + ,268670 + ,275960 + ,287040 + ,291500 + ,287540 + ,272340 + ,256570 + ,246190 + ,237340 + ,235890 + ,232790 + ,226100 + ,221110 + ,221950 + ,222990 + ,232840 + ,231380 + ,221690 + ,210190 + ,202210 + ,195210 + ,194050 + ,192360 + ,187170 + ,185930 + ,188840 + ,189350 + ,194320 + ,195850 + ,188840 + ,182380 + ,177200 + ,169140 + ,171900 + ,167320 + ,164130 + ,162550 + ,166260 + ,166010 + ,177070 + ,180800 + ,179570 + ,178130 + ,178180 + ,179430 + ,179130 + ,182220 + ,183650 + ,187110 + ,197780 + ,204000 + ,221330 + ,232500 + ,237570 + ,241100 + ,244420 + ,247240 + ,254300 + ,258970 + ,262230 + ,268880 + ,277950 + ,286630 + ,299830 + ,308090 + ,311400 + ,313520 + ,310980 + ,314910 + ,310950 + ,310720 + ,310380 + ,310570 + ,314790 + ,323400 + ,335010 + ,341630 + ,340820 + ,336280 + ,325490 + ,323750 + ,317510 + ,313890 + ,308610 + ,303720 + ,303090 + ,305140 + ,304040 + ,307100 + ,304330 + ,294710 + ,286890 + ,279050 + ,271860 + ,266710 + ,259590 + ,253830 + ,250640 + ,249140 + ,250840 + ,247590 + ,237830 + ,226380 + ,217230 + ,211420 + ,207620 + ,204310 + ,197490 + ,193580 + ,192330 + ,191970 + ,196070 + ,191940 + ,185620 + ,179410 + ,173920 + ,169190 + ,166840 + ,165170 + ,161450 + ,160830 + ,163670 + ,170830 + ,182690 + ,190940 + ,197770 + ,205090 + ,210720 + ,220210 + ,229730 + ,237070 + ,241620 + ,250370 + ,258570 + ,269860 + ,283220 + ,289610 + ,281770 + ,274700 + ,267650 + ,261380 + ,260500 + ,260730 + ,254200 + ,250450 + ,253380 + ,263740 + ,276240 + ,273820 + ,265890 + ,258400 + ,253520 + ,250710 + ,252850 + ,255260 + ,251170 + ,252500 + ,257780 + ,269900 + ,291590 + ,298870 + ,295570 + ,292100 + ,290870 + ,290580 + ,297970 + ,304010 + ,304340 + ,309850 + ,322320 + ,340170 + ,369280 + ,376690 + ,379700 + ,379520 + ,377770 + ,381560 + ,394580 + ,399320 + ,400370 + ,408200 + ,419070 + ,437730) > 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] 256950.1 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 4447.166 > (armose <- arm / armse) [1] 57.77839 > (geo <- geomean(x)) [1] 249707.3 > (har <- harmean(x)) [1] 242625.4 > (qua <- quamean(x)) [1] 264198.4 > (win <- winmean(x)) [,1] [,2] [1,] 256856.1 4427.134 [2,] 256754.3 4404.793 [3,] 256649.5 4381.322 [4,] 256637.2 4376.514 [5,] 256540.8 4352.744 [6,] 256160.2 4284.885 [7,] 256101.5 4273.552 [8,] 256118.2 4269.764 [9,] 256058.6 4254.955 [10,] 256097.2 4236.282 [11,] 255675.5 4174.080 [12,] 254882.4 4040.976 [13,] 254676.6 3997.537 [14,] 254234.0 3909.496 [15,] 254392.6 3873.053 [16,] 254335.9 3863.985 [17,] 254360.7 3848.681 [18,] 254000.7 3806.268 [19,] 253969.0 3782.335 [20,] 253006.5 3672.369 [21,] 252818.4 3652.631 [22,] 252794.3 3646.873 [23,] 252812.3 3618.388 [24,] 252526.0 3554.048 [25,] 252403.6 3538.018 [26,] 252093.5 3500.160 [27,] 252211.7 3484.613 [28,] 252367.7 3442.923 [29,] 252358.6 3432.955 [30,] 252211.8 3383.709 [31,] 252153.6 3376.535 [32,] 252427.0 3348.217 [33,] 252387.4 3344.604 [34,] 252451.2 3333.211 [35,] 252706.4 3301.486 [36,] 252794.5 3274.102 [37,] 252561.4 3251.855 [38,] 252529.7 3235.640 [39,] 252334.7 3217.044 [40,] 252180.5 3156.150 [41,] 252110.1 3131.437 [42,] 252166.9 3125.518 [43,] 252301.3 3100.578 [44,] 252441.1 3085.866 [45,] 252464.5 3078.458 [46,] 252764.0 3042.395 [47,] 252678.3 3022.380 [48,] 251865.8 2951.878 [49,] 252751.4 2823.754 [50,] 252983.2 2760.235 [51,] 252428.0 2698.683 [52,] 252406.4 2659.537 [53,] 252765.2 2566.608 [54,] 253099.9 2468.250 [55,] 253105.6 2442.374 [56,] 253283.5 2421.692 [57,] 254821.4 2252.631 [58,] 255634.0 2167.407 [59,] 255612.4 2118.681 [60,] 255034.3 2058.890 [61,] 254989.8 2035.980 [62,] 255025.4 2024.803 [63,] 255343.6 1993.818 [64,] 256317.0 1901.372 > (tri <- trimean(x)) [,1] [,2] [1,] 256504.5 4362.170 [2,] 256145.4 4292.693 [3,] 255831.1 4230.802 [4,] 255546.4 4173.591 [5,] 255258.7 4113.786 [6,] 254985.2 4055.578 [7,] 254774.0 4007.545 [8,] 254567.1 3958.186 [9,] 254353.2 3905.882 [10,] 254141.6 3852.031 [11,] 253920.8 3796.823 [12,] 253738.5 3745.899 [13,] 253628.2 3707.818 [14,] 253533.8 3671.904 [15,] 253474.5 3642.895 [16,] 253401.1 3615.328 [17,] 253330.1 3586.429 [18,] 253255.4 3556.640 [19,] 253203.8 3528.412 [20,] 253153.0 3499.988 [21,] 253162.3 3478.855 [22,] 253183.6 3457.449 [23,] 253206.8 3434.527 [24,] 253229.7 3411.819 [25,] 253269.4 3392.131 [26,] 253316.9 3371.726 [27,] 253382.3 3352.184 [28,] 253443.5 3331.824 [29,] 253498.6 3312.554 [30,] 253555.8 3291.952 [31,] 253621.9 3272.835 [32,] 253693.0 3252.086 [33,] 253753.3 3231.195 [34,] 253817.3 3208.184 [35,] 253880.6 3183.446 [36,] 253934.2 3158.388 [37,] 253985.8 3132.576 [38,] 254049.5 3105.449 [39,] 254116.8 3076.400 [40,] 254195.2 3045.353 [41,] 254283.1 3015.474 [42,] 254377.3 2983.910 [43,] 254472.6 2948.926 [44,] 254565.9 2911.734 [45,] 254656.8 2871.162 [46,] 254750.3 2826.071 [47,] 254834.9 2778.570 [48,] 254926.7 2726.672 [49,] 255056.9 2674.288 [50,] 255155.1 2627.720 [51,] 255247.8 2580.954 [52,] 255368.4 2533.356 [53,] 255495.6 2482.726 [54,] 255613.3 2434.520 [55,] 255722.3 2389.589 [56,] 255836.5 2340.277 [57,] 255948.7 2285.350 [58,] 255998.7 2241.924 [59,] 256015.0 2201.139 [60,] 256033.2 2158.877 [61,] 256078.9 2116.122 [62,] 256129.3 2068.156 [63,] 256181.1 2012.353 [64,] 256220.9 1950.029 > (midr <- midrange(x)) [1] 299280 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 254337.4 254926.7 254337.4 254926.7 254926.7 254337.4 254926.7 254834.9 > postscript(file="/var/wessaorg/rcomp/tmp/1fx7q1414050558.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/2f18z1414050558.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/3019b1414050558.tab") > > try(system("convert tmp/1fx7q1414050558.ps tmp/1fx7q1414050558.png",intern=TRUE)) character(0) > try(system("convert tmp/2f18z1414050558.ps tmp/2f18z1414050558.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.586 0.185 1.776