R version 2.12.2 (2011-02-25) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-redhat-linux-gnu (64-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(131.676 + ,135.050 + ,129.070 + ,137.792 + ,139.762 + ,142.917 + ,144.198 + ,142.648 + ,152.170 + ,136.022 + ,138.142 + ,138.135 + ,135.027 + ,132.911 + ,133.976 + ,137.012 + ,119.610 + ,118.106 + ,120.383 + ,133.185 + ,131.416 + ,134.248 + ,134.397 + ,127.728 + ,131.837 + ,125.955 + ,134.187 + ,143.291 + ,145.074 + ,149.812 + ,144.668 + ,147.253 + ,145.568 + ,155.564 + ,155.872 + ,156.323 + ,158.010 + ,155.598 + ,154.785 + ,157.294 + ,162.938 + ,157.283 + ,166.074 + ,169.282 + ,172.552 + ,174.055 + ,175.409 + ,173.696 + ,171.283 + ,173.322 + ,170.717 + ,174.229 + ,175.339 + ,173.511 + ,175.839 + ,173.816 + ,173.990 + ,174.777 + ,174.819 + ,176.726 + ,176.199 + ,180.952 + ,176.663 + ,182.346 + ,180.605 + ,182.497 + ,187.856 + ,190.020 + ,190.108 + ,193.288 + ,193.230 + ,199.068 + ,195.076 + ,191.563 + ,191.067 + ,186.665 + ,185.508 + ,184.371 + ,183.046 + ,175.714 + ,175.768 + ,171.029 + ,170.465 + ,170.102 + ,156.389 + ,124.291 + ,99.360 + ,86.675 + ,85.056 + ,128.236 + ,164.257 + ,162.401 + ,152.779 + ,156.005 + ,153.387 + ,153.190 + ,148.840 + ,144.211 + ,145.953 + ,145.542 + ,150.271 + ,147.489 + ,143.824 + ,134.754 + ,131.736 + ,126.304 + ,125.511 + ,125.495 + ,130.133 + ,126.257 + ,110.323 + ,98.417 + ,105.749 + ,120.665 + ,124.075 + ,127.245 + ,146.731 + ,144.979 + ,148.210 + ,144.670 + ,142.970 + ,142.524 + ,146.142 + ,146.522 + ,148.128 + ,148.798 + ,150.181 + ,152.388 + ,155.694 + ,160.662 + ,155.520 + ,158.262 + ,154.338 + ,158.196 + ,160.371 + ,154.856 + ,150.636 + ,145.899 + ,141.242 + ,140.834 + ,141.119 + ,139.104 + ,134.437 + ,129.425 + ,123.155 + ,119.273 + ,120.472 + ,121.523 + ,121.983 + ,123.658 + ,124.794 + ,124.827 + ,120.382 + ,117.395 + ,115.790 + ,114.283 + ,117.271 + ,117.448 + ,118.764 + ,120.550 + ,123.554 + ,125.412 + ,124.182 + ,119.828 + ,115.361 + ,114.226 + ,115.214 + ,115.864 + ,114.276 + ,113.469 + ,114.883 + ,114.172 + ,111.225 + ,112.149 + ,115.618 + ,118.002 + ,121.382 + ,120.663) > 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] 145.2148 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 1.815371 > (armose <- arm / armse) [1] 79.99182 > (geo <- geomean(x)) [1] 143.1923 > (har <- harmean(x)) [1] 141.1541 > (qua <- quamean(x)) [1] 147.2096 > (win <- winmean(x)) [,1] [,2] [1,] 145.2015 1.8100637 [2,] 145.3133 1.7851174 [3,] 145.3283 1.7826318 [4,] 145.4344 1.7573837 [5,] 145.5489 1.7397718 [6,] 145.5470 1.7315912 [7,] 145.5799 1.7270665 [8,] 145.5419 1.7067904 [9,] 145.5173 1.6946983 [10,] 145.4553 1.6855594 [11,] 145.3881 1.6759020 [12,] 145.2993 1.6642593 [13,] 145.3030 1.6545612 [14,] 145.3172 1.6503608 [15,] 145.2121 1.6343804 [16,] 145.2040 1.6281601 [17,] 144.8499 1.5826085 [18,] 144.8511 1.5811082 [19,] 144.9517 1.5601843 [20,] 144.9252 1.5542277 [21,] 144.9231 1.5526615 [22,] 144.9849 1.5451114 [23,] 144.9589 1.5393700 [24,] 145.0382 1.5296376 [25,] 145.0366 1.5145897 [26,] 145.0797 1.5091959 [27,] 145.0297 1.4968455 [28,] 145.0894 1.4855905 [29,] 145.0790 1.4844093 [30,] 145.0647 1.4797771 [31,] 145.0574 1.4762099 [32,] 145.0444 1.4706647 [33,] 145.0098 1.4668035 [34,] 144.9996 1.4380474 [35,] 144.7778 1.4086986 [36,] 144.8195 1.3945989 [37,] 144.9983 1.3654242 [38,] 145.0297 1.3520313 [39,] 144.9729 1.3415510 [40,] 144.8823 1.3137334 [41,] 144.1681 1.2356794 [42,] 143.7651 1.1909968 [43,] 143.5679 1.1491112 [44,] 143.4434 1.1357770 [45,] 143.1516 1.0812180 [46,] 143.0979 1.0723670 [47,] 142.5452 1.0224248 [48,] 142.6471 1.0096508 [49,] 142.6791 0.9974510 [50,] 142.4912 0.9788950 [51,] 142.7576 0.9536389 [52,] 142.6376 0.9183963 [53,] 142.7692 0.9029172 [54,] 142.9257 0.8718733 [55,] 142.9943 0.8585689 [56,] 143.1610 0.8340445 [57,] 143.5411 0.7955969 [58,] 143.6148 0.7873715 [59,] 143.6201 0.7844183 > (tri <- trimean(x)) [,1] [,2] [1,] 145.2507 1.7774820 [2,] 145.3010 1.7424793 [3,] 145.2945 1.7188470 [4,] 145.2828 1.6944457 [5,] 145.2426 1.6757272 [6,] 145.1769 1.6598584 [7,] 145.1100 1.6444687 [8,] 145.0362 1.6286971 [9,] 144.9659 1.6150300 [10,] 144.8969 1.6020712 [11,] 144.8331 1.5893213 [12,] 144.7748 1.5767679 [13,] 144.7237 1.5645869 [14,] 144.6708 1.5524712 [15,] 144.6152 1.5397481 [16,] 144.5667 1.5276253 [17,] 144.5175 1.5150437 [18,] 144.4930 1.5058780 [19,] 144.4677 1.4959459 [20,] 144.4348 1.4869572 [21,] 144.4027 1.4775820 [22,] 144.3698 1.4673715 [23,] 144.3321 1.4567568 [24,] 144.2948 1.4455386 [25,] 144.2517 1.4339649 [26,] 144.2074 1.4224491 [27,] 144.1592 1.4101307 [28,] 144.1122 1.3975334 [29,] 144.0604 1.3844788 [30,] 144.0074 1.3700641 [31,] 143.9533 1.3544336 [32,] 143.8977 1.3373467 [33,] 143.8408 1.3188029 [34,] 143.7835 1.2984762 [35,] 143.7245 1.2783915 [36,] 143.6740 1.2586858 [37,] 143.6195 1.2378627 [38,] 143.5545 1.2170942 [39,] 143.4854 1.1949332 [40,] 143.4161 1.1709466 [41,] 143.3481 1.1466131 [42,] 143.3103 1.1271448 [43,] 143.2893 1.1094964 [44,] 143.2765 1.0934843 [45,] 143.2688 1.0764583 [46,] 143.2742 1.0623576 [47,] 143.2742 1.0469259 [48,] 143.3164 1.0340035 [49,] 143.3474 1.0202448 [50,] 143.3785 1.0054244 [51,] 143.3785 0.9899516 [52,] 143.4514 0.9747020 [53,] 143.4900 0.9606319 [54,] 143.5246 0.9457291 [55,] 143.5246 0.9317367 [56,] 143.5811 0.9165412 [57,] 143.6020 0.9014158 [58,] 143.6050 0.8882273 [59,] 143.6050 0.8730882 > (midr <- midrange(x)) [1] 142.062 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 143.0616 143.2765 143.2765 143.2765 143.2688 143.2765 143.2765 143.2765 > postscript(file="/var/www/wessaorg/rcomp/tmp/16br71301583052.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/www/wessaorg/rcomp/tmp/2j1oc1301583052.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/www/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/www/wessaorg/rcomp/tmp/3g0aj1301583052.tab") > > try(system("convert tmp/16br71301583052.ps tmp/16br71301583052.png",intern=TRUE)) character(0) > try(system("convert tmp/2j1oc1301583052.ps tmp/2j1oc1301583052.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.270 0.140 1.431