R version 2.8.0 (2008-10-20) Copyright (C) 2008 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(205597 + ,205471 + ,211064 + ,212856 + ,217036 + ,219302 + ,219759 + ,221388 + ,220834 + ,221788 + ,222358 + ,222972 + ,224164 + ,224915 + ,226294 + ,224690 + ,227021 + ,229284 + ,229189 + ,230032 + ,229389 + ,231053 + ,232560 + ,232681 + ,231555 + ,231428 + ,232141 + ,234939 + ,235424 + ,235471 + ,236355 + ,238693 + ,236958 + ,237060 + ,239282 + ,238252 + ,241552 + ,236230 + ,238909 + ,240723 + ,242120 + ,242100 + ,243276 + ,244677 + ,243494 + ,244902 + ,245247 + ,245578 + ,243052 + ,238121 + ,241863 + ,241203 + ,243634 + ,242351 + ,245180 + ,246126 + ,244424 + ,245166 + ,247258 + ,245094 + ,246020 + ,243082 + ,245555 + ,243685 + ,247277 + ,245029 + ,246169 + ,246778 + ,244577 + ,246048 + ,245775 + ,245328 + ,245477 + ,241903 + ,243219 + ,248088 + ,248521 + ,247389 + ,249057 + ,248916 + ,249193 + ,250768 + ,253106 + ,249829 + ,249447 + ,246755 + ,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) > 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] 254924.8 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 1386.871 > (armose <- arm / armse) [1] 183.8130 > (geo <- geomean(x)) [1] 254195.1 > (har <- harmean(x)) [1] 253455.7 > (qua <- quamean(x)) [1] 255644.4 > (win <- winmean(x)) [,1] [,2] [1,] 254870.3 1377.3100 [2,] 254881.1 1360.5270 [3,] 254906.7 1355.5753 [4,] 254985.2 1340.8760 [5,] 255041.8 1332.0294 [6,] 255046.8 1328.7750 [7,] 255078.6 1322.3939 [8,] 255066.4 1314.5649 [9,] 255081.6 1311.5900 [10,] 255089.9 1304.8600 [11,] 255112.5 1298.6104 [12,] 255038.9 1270.1293 [13,] 254958.9 1251.5104 [14,] 254901.0 1240.6675 [15,] 254982.7 1224.2042 [16,] 254965.2 1207.8339 [17,] 255051.9 1173.2687 [18,] 255048.9 1170.9056 [19,] 254849.5 1147.0493 [20,] 254798.7 1127.1622 [21,] 254768.2 1100.2457 [22,] 254715.2 1086.1188 [23,] 254670.7 1078.7609 [24,] 254738.7 1070.1086 [25,] 254783.6 1063.1998 [26,] 254782.8 1059.7916 [27,] 255067.0 1022.9429 [28,] 255107.3 1012.8236 [29,] 255027.7 1004.0055 [30,] 255120.5 989.6472 [31,] 255121.8 985.8888 [32,] 255128.1 967.3186 [33,] 255047.2 956.7271 [34,] 255231.2 938.1583 [35,] 255250.3 935.4566 [36,] 255323.2 926.7626 [37,] 255340.0 920.6294 [38,] 255373.4 910.1479 [39,] 255652.9 882.4428 [40,] 255742.9 872.7499 [41,] 255791.0 863.8940 [42,] 255843.1 856.6116 [43,] 255827.4 853.6274 [44,] 255836.3 846.5304 [45,] 255815.2 843.8241 [46,] 255854.8 837.7237 [47,] 255965.2 817.9229 [48,] 255953.2 815.5644 [49,] 255906.2 805.4076 [50,] 255897.9 802.1366 [51,] 255920.2 794.2411 [52,] 255924.2 788.1608 [53,] 255815.2 776.2420 [54,] 255950.8 753.1230 [55,] 255964.2 746.9798 [56,] 255942.1 740.2047 [57,] 255874.4 723.3054 [58,] 255748.7 706.2737 [59,] 255665.7 695.9962 [60,] 255637.0 689.9053 [61,] 255563.0 682.9938 [62,] 255445.7 669.7550 [63,] 255410.3 662.5283 [64,] 255362.3 650.5667 > (tri <- trimean(x)) [,1] [,2] [1,] 254919.8 1351.0584 [2,] 254970.4 1322.9527 [3,] 255016.5 1302.2407 [4,] 255054.6 1282.0036 [5,] 255072.9 1264.7338 [6,] 255079.6 1248.4307 [7,] 255085.5 1231.6751 [8,] 255086.5 1214.9225 [9,] 255089.3 1198.2859 [10,] 255090.3 1180.9060 [11,] 255090.3 1163.2170 [12,] 255088.0 1145.0552 [13,] 255092.8 1129.1399 [14,] 255104.8 1114.2411 [15,] 255122.1 1099.4331 [16,] 255133.2 1085.2833 [17,] 255146.0 1071.7419 [18,] 255152.8 1060.5997 [19,] 255160.0 1048.7872 [20,] 255180.6 1038.2344 [21,] 255205.1 1028.5985 [22,] 255232.1 1020.4937 [23,] 255263.0 1012.8571 [24,] 255297.3 1005.1234 [25,] 255328.8 997.3984 [26,] 255358.7 989.5164 [27,] 255389.5 981.1569 [28,] 255406.3 974.9661 [29,] 255421.6 968.9197 [30,] 255441.4 962.8946 [31,] 255457.2 957.3200 [32,] 255473.4 951.3992 [33,] 255489.9 946.2052 [34,] 255510.6 941.1609 [35,] 255523.6 936.8907 [36,] 255536.1 932.2445 [37,] 255545.7 927.6298 [38,] 255554.9 922.8415 [39,] 255562.9 918.1790 [40,] 255559.0 914.9132 [41,] 255550.9 911.7927 [42,] 255540.5 908.7565 [43,] 255527.5 905.6813 [44,] 255514.6 902.2381 [45,] 255500.9 898.6914 [46,] 255487.4 894.6849 [47,] 255471.8 890.4181 [48,] 255450.8 886.9244 [49,] 255429.4 882.8903 [50,] 255409.1 878.8953 [51,] 255388.3 874.3452 [52,] 255365.5 869.5507 [53,] 255341.5 864.3159 [54,] 255321.1 859.1367 [55,] 255293.8 854.9398 [56,] 255264.5 850.2955 [57,] 255234.8 845.2064 [58,] 255206.4 840.6109 [59,] 255182.1 836.5943 [60,] 255160.3 832.5115 [61,] 255138.5 827.8760 [62,] 255118.9 822.6938 [63,] 255103.5 817.5923 [64,] 255088.9 811.8347 > (midr <- midrange(x)) [1] 255402 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 255323.0 255450.8 255323.0 255450.8 255450.8 255323.0 255450.8 255471.8 > postscript(file="/var/www/html/freestat/rcomp/tmp/1oxgd1229252388.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/freestat/rcomp/tmp/2j9dy1229252388.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/33kte1229252388.tab") > > system("convert tmp/1oxgd1229252388.ps tmp/1oxgd1229252388.png") > system("convert tmp/2j9dy1229252388.ps tmp/2j9dy1229252388.png") > > > proc.time() user system elapsed 1.328 0.523 1.441