R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: x86_64-pc-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(1.0137 + ,0.9834 + ,0.9643 + ,0.947 + ,0.906 + ,0.9492 + ,0.9397 + ,0.9041 + ,0.8721 + ,0.8552 + ,0.8564 + ,0.8973 + ,0.9383 + ,0.9217 + ,0.9095 + ,0.892 + ,0.8742 + ,0.8532 + ,0.8607 + ,0.9005 + ,0.9111 + ,0.9059 + ,0.8883 + ,0.8924 + ,0.8833 + ,0.87 + ,0.8758 + ,0.8858 + ,0.917 + ,0.9554 + ,0.9922 + ,0.9778 + ,0.9808 + ,0.9811 + ,1.0014 + ,1.0183 + ,1.0622 + ,1.0773 + ,1.0807 + ,1.0848 + ,1.1582 + ,1.1663 + ,1.1372 + ,1.1139 + ,1.1222 + ,1.1692 + ,1.1702 + ,1.2286 + ,1.2613 + ,1.2646 + ,1.2262 + ,1.1985 + ,1.2007 + ,1.2138 + ,1.2266 + ,1.2176 + ,1.2218 + ,1.249 + ,1.2991 + ,1.3408 + ,1.3119 + ,1.3014 + ,1.3201 + ,1.2938 + ,1.2694 + ,1.2165 + ,1.2037 + ,1.2292 + ,1.2256 + ,1.2015 + ,1.1786 + ,1.1856 + ,1.2103 + ,1.1938 + ,1.202 + ,1.2271 + ,1.277 + ,1.265 + ,1.2684 + ,1.2811 + ,1.2727 + ,1.2611 + ,1.2881 + ,1.3213 + ,1.2999 + ,1.3074 + ,1.3242 + ,1.3516 + ,1.3511 + ,1.3419 + ,1.3716 + ,1.3622 + ,1.3896 + ,1.4227 + ,1.4684 + ,1.457 + ,1.4718 + ,1.4748 + ,1.5527 + ,1.5751 + ,1.5557 + ,1.5553 + ,1.577 + ,1.4975 + ,1.4369 + ,1.3322 + ,1.2732 + ,1.3449 + ,1.3239 + ,1.2785 + ,1.305 + ,1.319 + ,1.365 + ,1.4016 + ,1.4088 + ,1.4268 + ,1.4562 + ,1.4816 + ,1.4914 + ,1.4614 + ,1.4272 + ,1.3686 + ,1.3569 + ,1.3406 + ,1.2565 + ,1.2209 + ,1.277 + ,1.2894 + ,1.3067 + ,1.3898 + ,1.3661 + ,1.322 + ,1.336 + ,1.3649 + ,1.3999 + ,1.4442 + ,1.4349 + ,1.4388 + ,1.4264 + ,1.4343 + ,1.377 + ,1.3706 + ,1.3556 + ,1.3179 + ,1.2905 + ,1.3224 + ,1.3201 + ,1.3162 + ,1.2789 + ,1.2526 + ,1.2288 + ,1.24 + ,1.2856 + ,1.2974 + ,1.2828 + ,1.3119 + ,1.3288 + ,1.3359 + ,1.2964 + ,1.3026 + ,1.2982 + ,1.3189 + ,1.308 + ,1.331 + ,1.3348 + ,1.3635 + ,1.3493 + ,1.3704 + ,1.361 + ,1.3658 + ,1.3823 + ,1.3812 + ,1.3732 + ,1.3592 + ,1.3539 + ,1.3316 + ,1.2901 + ,1.2673 + ,1.2472 + ,1.2331) > ylimmax = '' > ylimmin = '' > main = 'Robustness of Central Tendency' > #'GNU S' R Code compiled by R2WASP v. 1.2.291 () > #Author: root > #To cite this work: Wessa, P., (2012), Central Tendency (v1.0.4) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_centraltendency.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > # > 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] 1.235656 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 0.01364086 > (armose <- arm / armse) [1] 90.58494 > (geo <- geomean(x)) [1] 1.220903 > (har <- harmean(x)) [1] 1.204847 > (qua <- quamean(x)) [1] 1.249061 > (win <- winmean(x)) [,1] [,2] [1,] 1.235657 0.013637648 [2,] 1.235454 0.013606412 [3,] 1.235519 0.013594439 [4,] 1.235668 0.013555379 [5,] 1.234193 0.013361882 [6,] 1.234060 0.013329146 [7,] 1.233741 0.013279345 [8,] 1.233772 0.013198347 [9,] 1.233747 0.013164618 [10,] 1.233697 0.013125179 [11,] 1.233496 0.013049793 [12,] 1.233229 0.013017499 [13,] 1.233525 0.012960425 [14,] 1.232841 0.012836662 [15,] 1.232691 0.012752360 [16,] 1.232682 0.012714154 [17,] 1.232502 0.012695911 [18,] 1.232792 0.012640421 [19,] 1.232212 0.012550385 [20,] 1.232823 0.012453230 [21,] 1.233324 0.012371740 [22,] 1.234901 0.012050051 [23,] 1.233304 0.011874727 [24,] 1.233317 0.011662419 [25,] 1.233387 0.011601575 [26,] 1.232823 0.011364429 [27,] 1.234128 0.011181074 [28,] 1.235093 0.010812592 [29,] 1.235399 0.010734983 [30,] 1.234749 0.010675740 [31,] 1.234491 0.010574841 [32,] 1.235771 0.010347352 [33,] 1.237274 0.010113078 [34,] 1.239559 0.009809461 [35,] 1.240104 0.009668584 [36,] 1.248384 0.008531229 [37,] 1.251426 0.008151937 [38,] 1.251975 0.008053228 [39,] 1.252842 0.007946275 [40,] 1.258997 0.007173161 [41,] 1.260592 0.006936796 [42,] 1.263812 0.006528901 [43,] 1.268398 0.005961384 [44,] 1.269816 0.005710988 [45,] 1.270216 0.005610154 [46,] 1.270037 0.005547639 [47,] 1.271630 0.005278578 [48,] 1.273363 0.005085326 [49,] 1.275106 0.004830762 [50,] 1.275189 0.004604489 [51,] 1.274962 0.004475128 [52,] 1.274876 0.004427221 [53,] 1.274964 0.004408700 [54,] 1.274094 0.004247929 [55,] 1.276080 0.004060782 [56,] 1.276827 0.003934615 [57,] 1.276858 0.003790862 [58,] 1.277019 0.003743626 [59,] 1.277904 0.003632203 [60,] 1.277471 0.003546301 > (tri <- trimean(x)) [,1] [,2] [1,] 1.235887 0.013488832 [2,] 1.236123 0.013329581 [3,] 1.236468 0.013176274 [4,] 1.236799 0.013016509 [5,] 1.237099 0.012856690 [6,] 1.237721 0.012731758 [7,] 1.238383 0.012604059 [8,] 1.239111 0.012475520 [9,] 1.239852 0.012350309 [10,] 1.240616 0.012220210 [11,] 1.241404 0.012085146 [12,] 1.242233 0.011948930 [13,] 1.242233 0.011805298 [14,] 1.243984 0.011656560 [15,] 1.244939 0.011508984 [16,] 1.245932 0.011357785 [17,] 1.246953 0.011197260 [18,] 1.248015 0.011024030 [19,] 1.249087 0.010840864 [20,] 1.250229 0.010649344 [21,] 1.251364 0.010449875 [22,] 1.252501 0.010239616 [23,] 1.253576 0.010041985 [24,] 1.254778 0.009840249 [25,] 1.256016 0.009638393 [26,] 1.256016 0.009421191 [27,] 1.258633 0.009202967 [28,] 1.259951 0.008980047 [29,] 1.261261 0.008769090 [30,] 1.262598 0.008541767 [31,] 1.264014 0.008292287 [32,] 1.265492 0.008021753 [33,] 1.266959 0.007742441 [34,] 1.268404 0.007453845 [35,] 1.269793 0.007163234 [36,] 1.271206 0.006848524 [37,] 1.272283 0.006629668 [38,] 1.273259 0.006425948 [39,] 1.274247 0.006205671 [40,] 1.275235 0.005967317 [41,] 1.275981 0.005789760 [42,] 1.276684 0.005614868 [43,] 1.277271 0.005463283 [44,] 1.277675 0.005354392 [45,] 1.278032 0.005256572 [46,] 1.278388 0.005155072 [47,] 1.278767 0.005044770 [48,] 1.279093 0.004947955 [49,] 1.279093 0.004858100 [50,] 1.279550 0.004782347 [51,] 1.279751 0.004718179 [52,] 1.279751 0.004656358 [53,] 1.280212 0.004588255 [54,] 1.280460 0.004509491 [55,] 1.280763 0.004434289 [56,] 1.280988 0.004368941 [57,] 1.281191 0.004306190 [58,] 1.281405 0.004248041 [59,] 1.281624 0.004182741 [60,] 1.281813 0.004118298 > (midr <- midrange(x)) [1] 1.2151 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 1.276804 1.278032 1.276804 1.278032 1.278032 1.276804 1.278032 1.277675 > postscript(file="/var/wessaorg/rcomp/tmp/1tlaq1424974895.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/238na1424974895.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/3sxzm1424974895.tab") > > try(system("convert tmp/1tlaq1424974895.ps tmp/1tlaq1424974895.png",intern=TRUE)) character(0) > try(system("convert tmp/238na1424974895.ps tmp/238na1424974895.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.665 0.159 1.835