R version 2.9.0 (2009-04-17) Copyright (C) 2009 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(277 + ,260.6 + ,291.6 + ,275.4 + ,275.3 + ,231.7 + ,238.8 + ,274.2 + ,277.8 + ,299.1 + ,286.6 + ,232.3 + ,294.1 + ,267.5 + ,309.7 + ,280.7 + ,287.3 + ,235.7 + ,256.4 + ,289 + ,290.8 + ,321.9 + ,291.8 + ,241.4 + ,295.5 + ,258.2 + ,306.1 + ,281.5 + ,283.1 + ,237.4 + ,274.8 + ,299.3 + ,300.4 + ,340.9 + ,318.8 + ,265.7 + ,322.7 + ,281.6 + ,323.5 + ,312.6 + ,310.8 + ,262.8 + ,273.8 + ,320 + ,310.3 + ,342.2 + ,320.1 + ,265.6 + ,327 + ,300.7 + ,346.4 + ,317.3 + ,326.2 + ,270.7 + ,278.2 + ,324.6 + ,321.8 + ,343.5 + ,354 + ,278.2 + ,330.2 + ,307.3 + ,375.9 + ,335.3 + ,339.3 + ,280.3 + ,293.7 + ,341.2 + ,345.1 + ,368.7 + ,369.4 + ,288.4 + ,341 + ,319.1 + ,374.2 + ,344.5 + ,337.3 + ,281 + ,282.2 + ,321 + ,325.4 + ,366.3 + ,380.3 + ,300.7 + ,359.3 + ,327.6 + ,383.6 + ,352.4 + ,329.4 + ,294.5 + ,333.5 + ,334.3 + ,358 + ,396.1 + ,387 + ,307.2 + ,363.9 + ,344.7 + ,397.6 + ,376.8 + ,337.1 + ,299.3 + ,323.1 + ,329.1 + ,347 + ,462 + ,436.5 + ,360.4 + ,415.5 + ,382.1 + ,432.2 + ,424.3 + ,386.7 + ,354.5 + ,375.8 + ,368 + ,402.4 + ,426.5 + ,433.3 + ,338.5 + ,416.8 + ,381.1 + ,445.7 + ,412.4 + ,394 + ,348.2 + ,380.1 + ,373.7 + ,393.6 + ,434.2 + ,430.7 + ,344.5 + ,411.9 + ,370.5 + ,437.3 + ,411.3 + ,385.5 + ,341.3 + ,384.2 + ,373.2 + ,415.8 + ,448.6 + ,454.3 + ,350.3 + ,419.1 + ,398 + ,456.1 + ,430.1 + ,399.8 + ,362.7 + ,384.9 + ,385.3 + ,432.3 + ,468.9 + ,442.7 + ,370.2 + ,439.4 + ,393.9 + ,468.7 + ,438.8 + ,430.1 + ,366.3 + ,391 + ,380.9 + ,431.4 + ,465.4 + ,471.5 + ,387.5 + ,446.4 + ,421.5 + ,504.8 + ,492.1 + ,421.3 + ,396.7 + ,428 + ,421.9 + ,465.6 + ,525.8 + ,499.9 + ,435.3 + ,479.5 + ,473 + ,554.4 + ,489.6 + ,462.2 + ,420.3) > 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] 358.0833 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 4.924848 > (armose <- arm / armse) [1] 72.70952 > (geo <- geomean(x)) [1] 351.8634 > (har <- harmean(x)) [1] 345.7112 > (qua <- quamean(x)) [1] 364.2948 > (win <- winmean(x)) [,1] [,2] [1,] 357.9328 4.893571 [2,] 357.7435 4.849093 [3,] 357.6919 4.832620 [4,] 357.5543 4.802539 [5,] 357.5570 4.783177 [6,] 357.7151 4.675785 [7,] 357.5382 4.634197 [8,] 357.5769 4.613702 [9,] 357.5575 4.585055 [10,] 357.6973 4.567030 [11,] 357.5199 4.542551 [12,] 357.6231 4.528301 [13,] 357.6231 4.476018 [14,] 357.8414 4.449958 [15,] 357.3978 4.387621 [16,] 357.2946 4.363622 [17,] 356.8194 4.297563 [18,] 356.6161 4.272188 [19,] 356.7081 4.247389 [20,] 356.4715 4.202517 [21,] 356.1441 4.157112 [22,] 356.0731 4.149447 [23,] 356.1473 4.103392 [24,] 356.0957 4.087235 [25,] 355.9747 4.066135 [26,] 355.8909 4.043019 [27,] 355.7747 4.027940 [28,] 355.7145 4.003347 [29,] 355.8392 3.987852 [30,] 356.2747 3.919533 [31,] 356.2747 3.896280 [32,] 356.3608 3.867589 [33,] 356.4672 3.857501 [34,] 356.4124 3.787123 [35,] 356.2806 3.744318 [36,] 355.8935 3.697824 [37,] 355.7941 3.615103 [38,] 355.7941 3.599454 [39,] 355.8360 3.587556 [40,] 355.8360 3.546559 [41,] 356.3651 3.448444 [42,] 355.8909 3.393651 [43,] 355.6597 3.371290 [44,] 355.8489 3.341015 [45,] 355.1715 3.262732 [46,] 355.0478 3.251025 [47,] 356.2608 3.115547 [48,] 354.2478 2.876865 [49,] 353.5892 2.813129 [50,] 353.7505 2.713309 [51,] 353.8054 2.689235 [52,] 353.6935 2.654929 [53,] 354.0355 2.595425 [54,] 354.7903 2.426160 [55,] 355.2043 2.386762 [56,] 355.2043 2.371411 [57,] 354.6833 2.279535 [58,] 353.6231 2.184189 [59,] 353.7500 2.147292 [60,] 353.9113 2.117980 [61,] 353.5505 2.082503 [62,] 353.7505 2.055130 > (tri <- trimean(x)) [,1] [,2] [1,] 357.7033 4.813105 [2,] 357.4687 4.726933 [3,] 357.3267 4.659461 [4,] 357.1994 4.593562 [5,] 357.1057 4.531936 [6,] 357.0092 4.470650 [7,] 356.8820 4.427326 [8,] 356.7794 4.388415 [9,] 356.6690 4.349937 [10,] 356.5584 4.312819 [11,] 356.4293 4.275303 [12,] 356.3154 4.238052 [13,] 356.3154 4.199575 [14,] 356.0589 4.164062 [15,] 355.9071 4.128458 [16,] 355.7870 4.096550 [17,] 355.6717 4.064363 [18,] 355.5880 4.035907 [19,] 355.5162 4.007320 [20,] 355.4363 3.978404 [21,] 355.3694 3.950877 [22,] 355.3211 3.924754 [23,] 355.2757 3.896767 [24,] 355.2246 3.869953 [25,] 355.1750 3.841824 [26,] 355.1750 3.812638 [27,] 355.0894 3.782413 [28,] 355.0531 3.750355 [29,] 355.0188 3.717030 [30,] 354.9770 3.681516 [31,] 354.9121 3.647820 [32,] 354.8451 3.612406 [33,] 354.7717 3.575491 [34,] 354.6907 3.535344 [35,] 354.6095 3.496619 [36,] 354.5316 3.457106 [37,] 354.4688 3.417027 [38,] 354.4082 3.379237 [39,] 354.3454 3.338235 [40,] 354.2783 3.293286 [41,] 354.2087 3.246455 [42,] 354.1127 3.202388 [43,] 354.0340 3.157727 [44,] 353.9622 3.109339 [45,] 353.8792 3.057336 [46,] 353.8223 3.006026 [47,] 353.7685 2.948822 [48,] 353.6589 2.896808 [49,] 353.6330 2.861340 [50,] 353.6349 2.826578 [51,] 353.6349 2.796016 [52,] 353.6298 2.762674 [53,] 353.6187 2.727171 [54,] 353.6000 2.691779 [55,] 353.5461 2.667789 [56,] 353.4703 2.642856 [57,] 353.3903 2.614331 [58,] 353.3300 2.590373 [59,] 353.3162 2.571951 [60,] 353.2955 2.552875 [61,] 353.2656 2.532075 [62,] 353.2516 2.510049 > (midr <- midrange(x)) [1] 393.05 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 352.6394 353.2632 353.2632 353.2632 353.7685 353.2632 353.2632 353.2632 > postscript(file="/var/www/html/rcomp/tmp/1ixbd1264413385.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/22pfo1264413385.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/3qms11264413385.tab") > > try(system("convert tmp/1ixbd1264413385.ps tmp/1ixbd1264413385.png",intern=TRUE)) character(0) > try(system("convert tmp/22pfo1264413385.ps tmp/22pfo1264413385.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.988 0.342 1.134