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(1160.9 + ,1470.1 + ,1653 + ,1390 + ,1550.7 + ,1852.1 + ,1797.9 + ,1256.4 + ,1835.5 + ,1793.2 + ,1874.8 + ,1948.6 + ,1819.9 + ,1823.6 + ,2042.7 + ,1705.9 + ,1876.8 + ,1919.5 + ,1803.4 + ,1380.4 + ,1812.1 + ,1856.4 + ,1934.3 + ,1977.5 + ,1813.7 + ,1937.1 + ,2178.7 + ,1699.9 + ,1943.8 + ,2151 + ,1927.7 + ,1418.2 + ,1899.8 + ,1974.9 + ,1918.7 + ,1818.3 + ,1751 + ,1888 + ,1947.7 + ,1818.1 + ,2110.8 + ,2009 + ,2124.2 + ,1542.5 + ,2044.1 + ,2324.3 + ,2055.6 + ,2137.9 + ,2140.7 + ,2301.5 + ,2412.6 + ,2477 + ,2311.8 + ,2326.6 + ,2745.3 + ,1937.2 + ,2668.7 + ,2636.1 + ,2327.6 + ,2584.6 + ,2282.8 + ,2337 + ,2499.3 + ,2389.6 + ,2327.1 + ,2556.6 + ,2580 + ,1831.4 + ,2382.1 + ,2237.9 + ,2117.7 + ,2240.9 + ,1946.1 + ,2149.6 + ,2690 + ,2171.1 + ,2358.6 + ,2841.4 + ,3064.6 + ,2037.3 + ,2799.9 + ,2852.3 + ,2541.2 + ,2910 + ,2694.6 + ,3081.7 + ,3648.2 + ,2823.3 + ,3670.1 + ,3027.6 + ,3578.5 + ,2655.1 + ,3835.3 + ,3766 + ,3716 + ,3531.7 + ,3194 + ,3442.2 + ,3610 + ,3105.5 + ,3428.4 + ,3489.7 + ,3679 + ,2596.1 + ,3110.7 + ,3401.7 + ,3431.8 + ,3383.1 + ,3797.5 + ,3860.5 + ,4054.1 + ,4044.9 + ,4402.4 + ,4046.4 + ,4329.7 + ,3204.2 + ,4037.2 + ,4678.4 + ,4174.6 + ,4151.4 + ,3874.7 + ,3568 + ,3431 + ,3733.2 + ,3278.3 + ,3583.7 + ,4060.3 + ,2979 + ,4078.4 + ,4002.1 + ,3542.4 + ,3928.2 + ,3626 + ,3998.4 + ,4413.5 + ,3853.1 + ,3920.5 + ,4616.2 + ,4332.7 + ,3362.7 + ,3855.4 + ,4087.1 + ,3860.2 + ,4018.1 + ,3627.7 + ,3996 + ,4420.7 + ,4386.5 + ,4631.8 + ,4875.3 + ,4549.3 + ,3933 + ,4963.3 + ,4419.7 + ,4646.7 + ,5000.2 + ,4302.2 + ,4432.1 + ,5125.5 + ,4299.9 + ,5145.6 + ,4537.8 + ,4880.9 + ,4136.9 + ,4668.8 + ,4818.5 + ,4933.9 + ,4524.4 + ,4676.1 + ,4911 + ,5745 + ,4483.7 + ,4772.7 + ,5021.4 + ,5535.6 + ,4736 + ,5001.7 + ,5486.3 + ,4958.5 + ,4756.4 + ,4763.3 + ,4951.2 + ,4692.8 + ,5301.2 + ,4929.7 + ,5223.3 + ,5708.3 + ,4018) > 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] 3215.95 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 85.75549 > (armose <- arm / armse) [1] 37.50139 > (geo <- geomean(x)) [1] 2993.406 > (har <- harmean(x)) [1] 2773.801 > (qua <- quamean(x)) [1] 3423.088 > (win <- winmean(x)) [,1] [,2] [1,] 3216.263 85.66129 [2,] 3215.745 85.22798 [3,] 3215.111 85.09697 [4,] 3211.773 84.48697 [5,] 3211.081 84.06257 [6,] 3210.912 83.49879 [7,] 3210.469 83.37390 [8,] 3210.393 82.38962 [9,] 3211.695 82.05509 [10,] 3211.934 82.01439 [11,] 3212.414 81.50843 [12,] 3214.801 81.21804 [13,] 3214.621 81.12979 [14,] 3213.743 80.94465 [15,] 3214.102 80.84207 [16,] 3212.646 80.64952 [17,] 3210.322 80.30798 [18,] 3209.805 80.24663 [19,] 3204.227 79.60176 [20,] 3199.748 79.04167 [21,] 3199.569 78.84914 [22,] 3199.241 78.71914 [23,] 3198.777 78.26857 [24,] 3193.811 77.64372 [25,] 3194.343 77.22199 [26,] 3194.301 77.16405 [27,] 3194.861 76.90998 [28,] 3193.327 76.41453 [29,] 3193.944 75.91836 [30,] 3191.582 75.65548 [31,] 3181.903 74.43642 [32,] 3181.069 74.14338 [33,] 3179.209 73.86958 [34,] 3171.866 73.15661 [35,] 3163.488 72.13323 [36,] 3161.746 71.88853 [37,] 3161.864 71.84163 [38,] 3160.793 71.70799 [39,] 3163.946 71.00213 [40,] 3161.116 70.63788 [41,] 3156.253 68.94479 [42,] 3161.905 68.32403 [43,] 3156.850 67.64305 [44,] 3156.639 67.56547 [45,] 3129.400 64.66771 [46,] 3137.230 62.98973 [47,] 3135.330 62.52819 [48,] 3124.274 61.30498 [49,] 3125.578 60.80391 [50,] 3121.509 60.33724 [51,] 3122.241 59.98780 [52,] 3120.498 59.77755 [53,] 3125.742 59.25266 [54,] 3125.713 58.88147 [55,] 3137.445 56.94256 [56,] 3138.309 56.86465 [57,] 3146.191 55.40060 [58,] 3150.819 54.82690 [59,] 3153.298 54.49728 [60,] 3137.181 52.52129 [61,] 3136.370 52.33342 [62,] 3133.996 52.11457 > (tri <- trimean(x)) [,1] [,2] [1,] 3213.402 84.88104 [2,] 3210.478 84.04794 [3,] 3210.478 83.39529 [4,] 3205.198 82.74597 [5,] 3203.462 82.22681 [6,] 3203.462 81.76922 [7,] 3200.200 81.38947 [8,] 3198.597 81.00059 [9,] 3196.966 80.73432 [10,] 3195.135 80.48988 [11,] 3193.232 80.22593 [12,] 3193.232 79.99740 [13,] 3188.954 79.77691 [14,] 3186.634 79.54116 [15,] 3184.330 79.29855 [16,] 3181.938 79.03899 [17,] 3179.595 78.76956 [18,] 3177.359 78.50249 [19,] 3175.100 78.21101 [20,] 3173.153 77.94480 [21,] 3171.440 77.69590 [22,] 3169.692 77.43264 [23,] 3167.913 77.14813 [24,] 3167.913 76.86593 [25,] 3164.539 76.59943 [26,] 3162.891 76.33164 [27,] 3161.196 76.03351 [28,] 3159.420 75.71708 [29,] 3157.669 75.39888 [30,] 3155.832 75.07738 [31,] 3155.832 74.73431 [32,] 3152.692 74.44008 [33,] 3151.325 74.12638 [34,] 3150.002 73.78957 [35,] 3148.977 73.46054 [36,] 3148.305 73.16188 [37,] 3147.689 72.83605 [38,] 3147.046 72.46473 [39,] 3146.428 72.04965 [40,] 3145.646 71.62984 [41,] 3144.960 71.17773 [42,] 3144.463 70.79408 [43,] 3143.697 70.39821 [44,] 3143.122 69.99380 [45,] 3142.533 69.52953 [46,] 3143.104 69.22881 [47,] 3143.360 69.00480 [48,] 3143.360 68.76439 [49,] 3144.554 68.56646 [50,] 3145.382 68.35545 [51,] 3145.382 68.12503 [52,] 3147.487 67.86244 [53,] 3148.677 67.54958 [54,] 3149.694 67.20944 [55,] 3150.764 66.82159 [56,] 3151.363 66.53310 [57,] 3151.955 66.16917 [58,] 3152.219 65.86031 [59,] 3152.284 65.51905 [60,] 3152.237 65.11234 [61,] 3152.952 64.81207 [62,] 3152.952 64.42797 > (midr <- midrange(x)) [1] 3452.95 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 3132.491 3143.360 3132.491 3143.360 3143.360 3132.491 3143.360 3143.104 > postscript(file="/var/www/html/rcomp/tmp/1je4e1228231705.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/2h2yk1228231705.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/342kq1228231706.tab") > > system("convert tmp/1je4e1228231705.ps tmp/1je4e1228231705.png") > system("convert tmp/2h2yk1228231705.ps tmp/2h2yk1228231705.png") > > > proc.time() user system elapsed 0.991 0.354 1.125