R version 2.6.0 (2007-10-03) Copyright (C) 2007 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(622.755354196211 + ,419.557064448606 + ,337.667064448595 + ,587.237064448597 + ,-78.5029355514001 + ,568.837064448604 + ,111.996766841012 + ,336.956766841012 + ,-410.723233158992 + ,506.376766841009 + ,-207.193233158988 + ,342.366766841013 + ,395.671231421387 + ,256.272941673806 + ,221.182941673811 + ,-171.54705832619 + ,-252.987058326191 + ,-506.447058326191 + ,-1043.78735593379 + ,-941.627355933786 + ,-841.907355933784 + ,-600.707355933784 + ,-876.677355933786 + ,-265.917355933785 + ,-632.912891353406 + ,-85.1111811009849 + ,-339.801181100985 + ,-361.331181100985 + ,-139.271181100985 + ,302.068818899015 + ,-360.67147870858 + ,645.48852129142 + ,410.608521291422 + ,-300.891478708579 + ,1072.53852129142 + ,-225.101478708580 + ,109.802985871800 + ,706.304696124221 + ,382.814696124222 + ,358.684696124221 + ,1231.64469612422 + ,593.68469612422 + ,1026.44439851663 + ,641.904398516626 + ,79.9243985166273 + ,249.724398516626 + ,720.354398516625 + ,-9.28560148337375 + ,-389.081136902994 + ,-315.179426650573 + ,-800.369426650573 + ,-673.499426650573 + ,-920.739426650575 + ,-1487.59942665057 + ,370.863251817773 + ,198.223251817774 + ,-9.05674818222588 + ,870.343251817775 + ,572.773251817774 + ,-117.366748182226 + ,593.537716398154 + ,-31.160573349427 + ,87.7494266505743 + ,798.919426650575 + ,-181.820573349426 + ,199.819426650573 + ,314.47912904298 + ,58.83912904298 + ,-547.640870957020 + ,49.7591290429804 + ,-771.71087095702 + ,-701.150870957021 + ,-141.946406376641 + ,-1102.04469612422 + ,-506.434696124219 + ,209.935303875781 + ,-826.90469612422 + ,-268.264696124222 + ,-677.204993731815 + ,-493.544993731816 + ,-136.324993731815 + ,-191.124993731815 + ,-1176.29499373181 + ,438.965006268185 + ,-431.030529151435 + ,-146.428818899016 + ,86.0811811009855 + ,-81.7488188990137 + ,-65.288818899014 + ,493.651181100984 + ,-478.489116506611 + ,-296.029116506610 + ,53.1908834933911 + ,305.190883493390 + ,-424.979116506609 + ,84.5808834933917 + ,-889.51465192623 + ,-80.7129416738088 + ,432.597058326191 + ,-769.932941673809 + ,828.527058326191 + ,1028.86705832619 + ,638.626760718596 + ,307.086760718597 + ,1529.50676071860 + ,-202.193239281403 + ,800.236760718596 + ,299.296760718597 + ,-377.998774701024 + ,378.502935551396 + ,98.5129355513976 + ,103.282935551397 + ,405.342935551397 + ,-924.617064448604 + ,97.7426379438026 + ,-457.297362056198 + ,-127.577362056198 + ,-686.477362056198 + ,290.952637943803 + ,153.612637943802 + ,1140.71710252418) > 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: Wessa, P., (2007), Central Tendency (v1.0.2) 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 > #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] 6.796767e-14 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 51.40381 > (armose <- arm / armse) [1] 1.322230e-15 > (geo <- geomean(x)) [1] NaN Warning message: In log(x) : NaNs produced > (har <- harmean(x)) [1] -611.2891 > (qua <- quamean(x)) [1] 563.1005 > (win <- winmean(x)) [,1] [,2] [1,] 0.1110940 50.28515 [2,] -0.1645638 49.75523 [3,] -0.4105450 49.17902 [4,] 1.5229603 48.34667 [5,] 2.1257550 48.21552 [6,] -5.4225181 46.89584 [7,] -6.0352448 46.23487 [8,] -7.0569309 45.82218 [9,] -4.5687161 45.40451 [10,] -9.8218044 44.28873 [11,] -8.6867528 43.74665 [12,] -11.8759382 42.52112 [13,] -12.0699920 42.44316 [14,] -4.4909667 41.24149 [15,] -4.6394664 40.73245 [16,] -7.2574221 40.07459 [17,] -6.7574553 39.99959 [18,] -1.6570759 39.05600 [19,] 1.1288144 38.09778 [20,] 9.2495248 36.88647 [21,] 5.5586472 34.65742 [22,] 3.2471520 34.38031 [23,] -4.6976320 32.87409 [24,] -2.9744064 32.36437 [25,] -1.2901583 31.53150 [26,] 2.4311271 30.64500 [27,] 2.6064769 30.35601 [28,] 3.6672787 29.72403 [29,] 5.7729090 28.78874 [30,] 7.4515706 28.35948 [31,] 9.7645062 27.66580 [32,] 6.7181978 27.29807 [33,] 7.9597528 26.18794 [34,] 13.5576848 25.28949 [35,] 17.4851043 24.82276 [36,] 12.2441959 23.95761 [37,] 18.4736664 22.81974 [38,] 18.6154481 22.67989 [39,] 21.7767794 22.14044 [40,] 30.0787700 21.10196 > (tri <- trimean(x)) [,1] [,2] [1,] -0.3521625 49.07243 [2,] -0.8312567 47.72600 [3,] -1.1819952 46.54180 [4,] -1.4573505 45.46273 [5,] -2.2695523 44.53106 [6,] -3.2453912 43.52730 [7,] -2.8350605 42.71942 [8,] -2.3082274 41.95451 [9,] -1.6109056 41.17599 [10,] -1.2171817 40.37640 [11,] -0.1655056 39.66338 [12,] 0.8008214 38.94483 [13,] 2.1463371 38.31697 [14,] 3.5691459 37.61433 [15,] 4.3346668 36.98888 [16,] 5.1480527 36.34482 [17,] 6.2264021 35.69594 [18,] 7.3136317 34.95997 [19,] 8.0401750 34.24227 [20,] 8.5835633 33.54272 [21,] 8.5325624 32.89422 [22,] 8.7551003 32.42956 [23,] 9.1590165 31.91728 [24,] 10.1576195 31.50324 [25,] 11.0901167 31.07415 [26,] 11.9585302 30.66735 [27,] 12.6203073 30.29170 [28,] 13.3107196 29.87078 [29,] 13.9722027 29.44174 [30,] 14.5330363 29.04334 [31,] 14.5330363 28.60837 [32,] 15.3768242 28.16676 [33,] 15.9721048 27.66523 [34,] 16.5264185 27.20674 [35,] 16.7335792 26.76585 [36,] 16.6805562 26.27221 [37,] 16.9978136 25.78053 [38,] 16.8905595 25.34027 [39,] 16.7628290 24.77270 [40,] 16.3834119 24.12047 > (midr <- midrange(x)) [1] 20.95367 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 8.466871 14.533036 14.533036 14.533036 14.533036 8.023130 14.533036 [8] 14.533036 > postscript(file="/var/www/html/rcomp/tmp/1hexf1198253050.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/2gisg1198253050.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 > 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/32wy81198253051.tab") > > system("convert tmp/1hexf1198253050.ps tmp/1hexf1198253050.png") > system("convert tmp/2gisg1198253050.ps tmp/2gisg1198253050.png") > > > proc.time() user system elapsed 1.118 0.329 1.223