R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-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(186.59 + ,244.665 + ,248.18 + ,253.568 + ,171.242 + ,413.971 + ,216.89 + ,227.901 + ,259.823 + ,148.438 + ,241.013 + ,206.248 + ,108.908 + ,267.952 + ,314.219 + ,235.115 + ,203.027 + ,365.415 + ,350.933 + ,263.304 + ,738.751 + ,959.073 + ,483.828 + ,213.016 + ,177.341 + ,352.622 + ,352.622 + ,217.307 + ,236.184 + ,215.701 + ,228.383 + ,485.625 + ,252.502 + ,342.515 + ,196.931 + ,365.315 + ,316.664 + ,313.523 + ,188.124 + ,184.083 + ,362.962 + ,170.161 + ,167.484 + ,211.752 + ,276.469 + ,182.097 + ,266.904 + ,235.328 + ,329.45 + ,258.118 + ,952.515 + ,247.141 + ,498.726 + ,313.308 + ,420.188 + ,231.156 + ,227.844 + ,204.385 + ,216.563 + ,207.19 + ,232.417 + ,203.109 + ,221.07 + ,254.623 + ,108.981 + ,229.417 + ,476.376 + ,221.91 + ,158.946 + ,295.746 + ,187.976 + ,283.76 + ,705.401 + ,178.69 + ,232.635 + ,245.185 + ,186.03 + ,181.142 + ,228.478 + ,491.849 + ,506.461 + ,182.828 + ,263.654 + ,253.718 + ,480.937 + ,209.894 + ,655.92 + ,223.408 + ,103.138 + ,255.664 + ,184.661 + ,372.945 + ,758.6 + ,256.445 + ,204.868 + ,197.384 + ,245.311 + ,301.707 + ,501.478 + ,278.731 + ,205.92 + ,177.14 + ,139.753 + ,366.46 + ,435.522 + ,239.906 + ,178.722 + ,340.04 + ,236.948 + ,221.152 + ,263.303 + ,222.814 + ,317.99 + ,149.593 + ,221.983 + ,201.551 + ,266.901 + ,185.218 + ,347.536 + ,395.593 + ,238.217 + ,254.697 + ,157.584 + ,807.302 + ,252.391 + ,189.194 + ,267.834 + ,173.15 + ,267.633 + ,283.284 + ,209.475 + ,135.135 + ,285.012 + ,178.495 + ,256.852 + ,301.828 + ,158.403 + ,355.963 + ,364.117 + ,233.203 + ,257.634 + ,208.57 + ,212.503 + ,251.059 + ,276.803 + ,198.339 + ,301.018 + ,369.761 + ,162.768 + ,199.968 + ,406.676 + ,364.156 + ,202.391 + ,319.491 + ,185.546 + ,243.559 + ,220.928 + ,193.714 + ,372.943 + ,163.822 + ,217.566 + ,232.527) > 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] 281.0759 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 11.02478 > (armose <- arm / armse) [1] 25.49493 > (geo <- geomean(x)) [1] 257.7475 > (har <- harmean(x)) [1] 240.7068 > (qua <- quamean(x)) [1] 313.9626 > (win <- winmean(x)) [,1] [,2] [1,] 281.0711 11.005859 [2,] 279.2792 10.381273 [3,] 278.8617 10.058870 [4,] 278.4856 9.905522 [5,] 277.7243 9.591468 [6,] 275.9345 9.095872 [7,] 269.8217 7.582142 [8,] 269.6160 7.531048 [9,] 269.4933 7.499151 [10,] 269.3047 7.397835 [11,] 268.9537 7.313067 [12,] 269.0918 7.264706 [13,] 269.0747 7.203727 [14,] 268.7739 7.124098 [15,] 265.1678 6.454804 [16,] 264.0474 6.177000 [17,] 263.4161 6.073704 [18,] 262.7338 5.938989 [19,] 261.4568 5.744803 [20,] 258.6647 5.355798 [21,] 258.9781 5.326986 [22,] 258.6757 5.258171 [23,] 258.3108 5.187693 [24,] 258.3419 5.150950 [25,] 258.4157 5.140989 [26,] 258.3191 5.109055 [27,] 258.3672 5.103367 [28,] 258.2513 5.070309 [29,] 257.0986 4.902544 [30,] 256.7366 4.802552 [31,] 256.7649 4.800035 [32,] 256.6426 4.740006 [33,] 256.8714 4.574472 [34,] 256.4928 4.388017 [35,] 256.0559 4.314879 [36,] 253.9148 4.016951 [37,] 252.0123 3.723804 [38,] 252.0315 3.652184 [39,] 251.9145 3.599122 [40,] 251.4678 3.518527 [41,] 251.3124 3.497264 [42,] 251.5875 3.462934 [43,] 248.6686 3.121727 [44,] 248.9214 3.093462 [45,] 248.8211 3.065335 [46,] 247.5916 2.885369 [47,] 244.8778 2.535124 [48,] 244.7750 2.474932 [49,] 244.7578 2.449641 [50,] 243.9260 2.263164 [51,] 243.5554 2.184607 [52,] 243.6129 2.159962 [53,] 241.7049 1.825848 [54,] 241.9529 1.797015 > (tri <- trimean(x)) [,1] [,2] [1,] 277.9506 10.258548 [2,] 274.7511 9.400347 [3,] 272.3999 8.817765 [4,] 270.1341 8.302194 [5,] 267.9089 7.772837 [6,] 265.7887 7.267173 [7,] 263.9378 6.827440 [8,] 263.0051 6.677245 [9,] 262.0755 6.522114 [10,] 261.1352 6.357232 [11,] 260.1898 6.192141 [12,] 259.2546 6.023327 [13,] 258.2781 5.843520 [14,] 257.2740 5.652653 [15,] 256.2659 5.451593 [16,] 255.5264 5.320967 [17,] 254.8524 5.211377 [18,] 254.2047 5.102804 [19,] 253.5857 4.998565 [20,] 253.0356 4.905756 [21,] 252.6556 4.846023 [22,] 252.2423 4.782290 [23,] 251.8339 4.718436 [24,] 251.4337 4.654477 [25,] 251.0173 4.586597 [26,] 250.5815 4.511522 [27,] 250.1351 4.430625 [28,] 249.6691 4.340233 [29,] 249.1917 4.241956 [30,] 248.7587 4.149876 [31,] 248.3279 4.056606 [32,] 247.8780 3.951011 [33,] 247.4158 3.837555 [34,] 246.9219 3.726183 [35,] 246.4263 3.620132 [36,] 245.9310 3.507212 [37,] 245.5228 3.414736 [38,] 245.1924 3.343456 [39,] 244.8453 3.268930 [40,] 244.4872 3.188623 [41,] 244.1338 3.105218 [42,] 243.7701 3.010461 [43,] 243.3734 2.902901 [44,] 243.1038 2.825545 [45,] 242.8063 2.737349 [46,] 242.4970 2.636455 [47,] 242.2331 2.543996 [48,] 242.0950 2.486399 [49,] 241.9537 2.425489 [50,] 241.8042 2.355436 [51,] 241.6896 2.299307 [52,] 241.5874 2.242815 [53,] 241.4747 2.176926 [54,] 241.4617 2.150098 > (midr <- midrange(x)) [1] 531.1055 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 243.6263 244.4872 244.4872 244.4872 244.1338 244.4872 244.4872 244.4872 > postscript(file="/var/wessaorg/rcomp/tmp/1z7l11324145161.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/2bhrc1324145161.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/35uzm1324145161.tab") > > try(system("convert tmp/1z7l11324145161.ps tmp/1z7l11324145161.png",intern=TRUE)) character(0) > try(system("convert tmp/2bhrc1324145161.ps tmp/2bhrc1324145161.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.084 0.147 1.230