R version 2.7.0 (2008-04-22) 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(283.9 + ,340.5 + ,383 + ,335.4 + ,356 + ,431 + ,441.9 + ,262.3 + ,481.1 + ,442.8 + ,458.4 + ,472.7 + ,498.8 + ,463.3 + ,538.3 + ,494.2 + ,446.2 + ,497.5 + ,433.2 + ,359.4 + ,455.9 + ,483 + ,469.7 + ,423.8 + ,371 + ,435.7 + ,438.1 + ,364.5 + ,392.7 + ,405.1 + ,391.4 + ,282.8 + ,408 + ,480.9 + ,416.9 + ,350.7 + ,399.9 + ,438.1 + ,454.9 + ,407.7 + ,485 + ,494 + ,521.2 + ,333 + ,520.3 + ,570.5 + ,501.8 + ,408.4 + ,512.7 + ,596.9 + ,581.5 + ,605.7 + ,581.1 + ,538.3 + ,731.4 + ,457.7 + ,706.3 + ,708.5 + ,683.5 + ,722.1 + ,727.7 + ,674.4 + ,675.3 + ,628.7 + ,658.3 + ,786 + ,736.9 + ,567.3 + ,736.3 + ,719.6 + ,676.8 + ,589.5 + ,596.5 + ,661 + ,840.9 + ,683.1 + ,677.4 + ,780.2 + ,928.8 + ,557.5 + ,760.6 + ,847 + ,719 + ,699.8 + ,745.4 + ,871.5 + ,1199.8 + ,795.4 + ,1054.4 + ,753.7 + ,1093.8 + ,768.4 + ,1405.7 + ,1149.2 + ,1086.8 + ,994.4 + ,958.8 + ,1096.3 + ,1107.4 + ,862.4 + ,983 + ,1021.1 + ,1036.3 + ,684.5 + ,874.1 + ,1031.6 + ,1126.3 + ,1041.9 + ,1419.4 + ,1402.2 + ,1503.6 + ,1503.6 + ,1571.2 + ,1351.2 + ,1448.3 + ,1292.9 + ,1293.7 + ,1824.7 + ,1641.4 + ,1529.1 + ,1417 + ,1060 + ,1117 + ,1226.2 + ,1124 + ,1246 + ,1398.1 + ,1095.2 + ,1507.1 + ,1467 + ,1262.1 + ,1262.9 + ,1184.9 + ,1412.9 + ,1581.8 + ,1348.8 + ,1443.7 + ,1628.7 + ,1393.2 + ,1108.1 + ,1074.8 + ,1323.5 + ,1367 + ,1190.7 + ,1149.2 + ,1281.9 + ,1406.8 + ,1392.4 + ,1674.5 + ,1529.8 + ,1180.9 + ,1282.8 + ,1524.2 + ,1524.6 + ,1520.1 + ,1881.9 + ,1375.2 + ,1336.5 + ,1653.9 + ,1404.3 + ,1731.4 + ,1330.8 + ,1616.3 + ,1503.7 + ,1476.7 + ,1565.6 + ,1581.8 + ,1407.7 + ,1393.7 + ,1662.4 + ,1929 + ,1352.8 + ,1339.8 + ,1251.2 + ,1491.7 + ,1540.7 + ,1383.7 + ,1827.5 + ,1591.3 + ,1289.2 + ,1291.9 + ,1368.3 + ,1264.6 + ,1448.8 + ,1239.2 + ,1280.9 + ,1269.4 + ,1068.5) > 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] 964.1543 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 32.27423 > (armose <- arm / armse) [1] 29.87381 > (geo <- geomean(x)) [1] 853.3758 > (har <- harmean(x)) [1] 744.1272 > (qua <- quamean(x)) [1] 1060.367 > (win <- winmean(x)) [,1] [,2] [1,] 964.0128 32.22262 [2,] 963.4457 32.13561 [3,] 964.1846 32.04362 [4,] 962.2505 31.76712 [5,] 960.8729 31.56336 [6,] 960.8122 31.48303 [7,] 960.6931 31.42513 [8,] 960.3059 31.34800 [9,] 959.9420 31.25294 [10,] 959.6282 31.14293 [11,] 958.8676 30.90997 [12,] 958.7973 30.79086 [13,] 958.8872 30.78201 [14,] 958.6340 30.64460 [15,] 958.6021 30.55667 [16,] 956.7043 30.31350 [17,] 955.7457 30.21002 [18,] 955.7170 30.19950 [19,] 956.1213 30.07061 [20,] 956.8128 29.99624 [21,] 957.1590 29.87391 [22,] 955.8952 29.69752 [23,] 955.7851 29.62756 [24,] 956.0787 29.59758 [25,] 956.0787 29.59758 [26,] 954.9585 29.38668 [27,] 952.9335 29.16595 [28,] 951.9952 28.98062 [29,] 950.5298 28.59156 [30,] 950.6096 28.56934 [31,] 950.1479 28.47139 [32,] 946.1309 28.08221 [33,] 946.5697 27.96459 [34,] 946.9856 27.79137 [35,] 946.5761 27.65358 [36,] 947.9739 27.49493 [37,] 947.7968 27.47205 [38,] 947.8979 27.41198 [39,] 947.8771 27.33563 [40,] 948.9197 27.08549 [41,] 948.0037 26.99666 [42,] 948.6293 26.92063 [43,] 948.7436 26.87784 [44,] 947.4096 26.63596 [45,] 947.9840 26.22573 [46,] 948.1553 25.91439 [47,] 948.0553 25.86637 [48,] 948.7957 25.17069 [49,] 948.3787 25.13492 [50,] 952.8468 24.63827 [51,] 953.0638 24.20167 [52,] 953.0362 24.04838 [53,] 954.4176 23.65833 [54,] 952.4356 23.47080 [55,] 946.0580 22.54477 [56,] 947.9048 22.34927 [57,] 947.7229 22.31406 [58,] 949.6048 22.01784 [59,] 954.8144 21.25447 [60,] 963.9739 20.46568 [61,] 964.5255 20.36883 [62,] 965.1521 19.70428 > (tri <- trimean(x)) [,1] [,2] [1,] 962.7403 31.98168 [2,] 961.4402 31.72405 [3,] 961.4402 31.49615 [4,] 959.0883 31.28618 [5,] 958.2534 31.14042 [6,] 958.2534 31.03072 [7,] 957.1322 30.92724 [8,] 956.5762 30.82432 [9,] 956.0606 30.72385 [10,] 955.5780 30.62730 [11,] 955.1193 30.53561 [12,] 955.1193 30.46292 [13,] 954.3352 30.39526 [14,] 953.9238 30.32013 [15,] 953.5234 30.25015 [16,] 953.1154 30.17993 [17,] 952.8416 30.12373 [18,] 952.6303 30.06858 [19,] 952.4153 30.00584 [20,] 952.1676 29.94499 [21,] 951.8685 29.88108 [22,] 951.5396 29.81753 [23,] 951.2775 29.75838 [24,] 951.2775 29.69501 [25,] 950.7268 29.62374 [26,] 950.4309 29.54152 [27,] 950.1866 29.46347 [28,] 950.0417 29.39020 [29,] 949.9408 29.31889 [30,] 949.9109 29.26419 [31,] 949.9109 29.20000 [32,] 949.8629 29.13089 [33,] 950.0426 29.07714 [34,] 950.2075 29.01995 [35,] 950.3585 28.96275 [36,] 950.5336 28.90261 [37,] 950.6509 28.84058 [38,] 950.7804 28.76622 [39,] 950.9100 28.68093 [40,] 951.0454 28.58463 [41,] 951.1396 28.48912 [42,] 951.2779 28.38177 [43,] 951.3941 28.26019 [44,] 951.5100 28.11995 [45,] 951.6888 27.97465 [46,] 951.8500 27.83712 [47,] 952.0106 27.69963 [48,] 952.0106 27.53977 [49,] 952.3300 27.40973 [50,] 952.5023 27.25546 [51,] 952.5023 27.11444 [52,] 952.4619 26.98217 [53,] 952.4366 26.83354 [54,] 952.3487 26.68809 [55,] 952.3449 26.52654 [56,] 952.6276 26.41763 [57,] 952.8419 26.29597 [58,] 953.0764 26.14253 [59,] 953.2371 25.97998 [60,] 953.1632 25.85720 [61,] 952.6500 25.77963 [62,] 952.6500 25.67665 > (midr <- midrange(x)) [1] 1095.65 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 947.4663 952.0106 947.4663 952.0106 952.0106 947.4663 952.0106 951.8500 > postscript(file="/var/www/html/rcomp/tmp/1cc0u1228233790.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/rcomp/tmp/23ih31228233790.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/rcomp/tmp/33j0g1228233791.tab") > > system("convert tmp/1cc0u1228233790.ps tmp/1cc0u1228233790.png") > system("convert tmp/23ih31228233790.ps tmp/23ih31228233790.png") > > > proc.time() user system elapsed 2.603 0.590 2.695