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(-0.0188617052236618 + ,-0.00480518696381852 + ,0.0386604299596439 + ,-0.108558199150346 + ,-0.0342341157949376 + ,0.0470709618031273 + ,-0.00378593260964411 + ,0.084908753384194 + ,0.0216717154148097 + ,0.0764535262723437 + ,-0.0214701349220553 + ,-0.00288912125847166 + ,-0.0110750040363763 + ,-0.0131796528909636 + ,-0.0252908488260806 + ,-0.0723209465011967 + ,0.035583254480575 + ,0.0682703692876727 + ,-0.0212556897253976 + ,-0.0449837816224749 + ,0.0387402666693091 + ,0.0505999046218331 + ,0.0124994452158954 + ,-0.0164536913959678 + ,0.00337538175968105 + ,-0.0310346549045335 + ,0.00144469116877099 + ,-0.0512908613724161 + ,-0.003825501183449 + ,0.0561771288521043 + ,-0.0361660789107442 + ,0.127052534137075 + ,-0.0430000497266376 + ,-0.00300482147163695 + ,-0.0284425492263112 + ,0.028552270866014 + ,0.0672346139508214 + ,-0.0231487087966499 + ,0.00160225316569137 + ,-0.0359537994734657 + ,0.0626134330196644 + ,-0.02262721585705 + ,-0.0281902975621727 + ,-0.0372562976633048 + ,0.0801283787641748 + ,0.0392978596059758 + ,0.0433743926024704 + ,0.00694555229977395 + ,0.00178164632939194 + ,-0.0556942121591666 + ,0.00236425589354194 + ,-0.0269163708685144 + ,0.00628712124442563 + ,0.0338350851264338 + ,0.0903095429199688 + ,-0.0104845911526373 + ,0.00309375706670886 + ,-0.0567838417208723 + ,0.0261608296190736 + ,0.00705041100420995 + ,0.0149776155524088 + ,-0.0163278603742414 + ,-0.0120185263685339 + ,-0.00638048949476752 + ,0.0385711103907229 + ,0.0732916611361229 + ,0.0918216918380306 + ,0.0391920798231936 + ,0.0792929171874763 + ,-0.0412888196720461 + ,-0.0126660759823600 + ,-0.0389329639109142 + ,-0.0391534633552554 + ,0.00865927349803817 + ,-0.0366550619725453 + ,0.00948017353160925 + ,0.0282403016997684 + ,-0.0245998776547136 + ,-0.0405062485787519 + ,0.0579704575275294 + ,-0.0212434550014502 + ,-0.00913346491345895 + ,0.0713951189633145 + ,0.0333026486107708 + ,0.0273541497799685 + ,0.00900169997204778 + ,0.046362769940699 + ,0.0214786912951184 + ,-0.0391502769651846 + ,0.00583943395571575 + ,-0.0689344693008103 + ,-0.045789995214643 + ,0.0199823053845984 + ,-0.0290597654043378 + ,0.00468433591433152 + ,0.00168585008920773 + ,-0.0321491302845718 + ,0.00876022617539829 + ,0.0321908668323309 + ,-0.0532497260755194 + ,0.036189226577314 + ,0.0244785447090322 + ,0.0639802307638794 + ,-0.0734871807180395 + ,0.0278000207980543 + ,-0.00819467423881322 + ,-0.0115617093712317 + ,0.0493591368131129 + ,0.0177729429363333 + ,0.05212603893352 + ,-0.0215389488063387 + ,0.0551258163993368 + ,0.0313657098241115 + ,-0.0162011713075296 + ,0.0360045363634907 + ,-0.0915895507553786 + ,0.0592615564207261 + ,0.105372426104259 + ,0.0488645877985787 + ,-0.00816410416000369 + ,-0.00955522793648373 + ,0.00879456045456353 + ,-0.0153263431374198 + ,-0.0662904810276389 + ,0.0237068904542412 + ,0.0886587441320772 + ,-0.0451837643288755 + ,0.0591687634000339 + ,0.0876070888157847) > 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] 0.007248025 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 0.003878686 > (armose <- arm / armse) [1] 1.868681 > (geo <- geomean(x)) [1] NaN Warning message: In log(x) : NaNs produced > (har <- harmean(x)) [1] 0.05430973 > (qua <- quamean(x)) [1] 0.04447687 > (win <- winmean(x)) [,1] [,2] [1,] 0.007211502 0.003812804 [2,] 0.007282070 0.003720736 [3,] 0.007274026 0.003709967 [4,] 0.007327845 0.003683784 [5,] 0.007389564 0.003660433 [6,] 0.007706229 0.003573111 [7,] 0.007505956 0.003522003 [8,] 0.007605741 0.003492731 [9,] 0.007544308 0.003442933 [10,] 0.007725626 0.003349924 [11,] 0.007615600 0.003318981 [12,] 0.007343528 0.003273707 [13,] 0.007439061 0.003233916 [14,] 0.007271587 0.003161188 [15,] 0.007203654 0.003128171 [16,] 0.006955704 0.003051659 [17,] 0.006943896 0.003049972 [18,] 0.006807013 0.003024156 [19,] 0.006789830 0.002960542 [20,] 0.006720051 0.002928586 [21,] 0.006311317 0.002857229 [22,] 0.006087249 0.002820671 [23,] 0.006172637 0.002758101 [24,] 0.006468532 0.002702983 [25,] 0.006336914 0.002636692 [26,] 0.006592217 0.002575906 [27,] 0.006095927 0.002487559 [28,] 0.005265851 0.002380318 [29,] 0.005528457 0.002346464 [30,] 0.005801413 0.002294248 [31,] 0.005948275 0.002274592 [32,] 0.006286098 0.002234687 [33,] 0.005810184 0.002152954 [34,] 0.006048336 0.002118188 [35,] 0.005952705 0.002103684 [36,] 0.005524689 0.002044336 [37,] 0.005375484 0.002027498 [38,] 0.005749583 0.001921155 [39,] 0.006228121 0.001822452 [40,] 0.005394754 0.001725873 [41,] 0.005335867 0.001711534 [42,] 0.005477347 0.001668677 [43,] 0.006044287 0.001584314 > (tri <- trimean(x)) [,1] [,2] [1,] 0.007216543 0.003713304 [2,] 0.007221745 0.003603299 [3,] 0.007190110 0.003535554 [4,] 0.007160289 0.003465157 [5,] 0.007114880 0.003395683 [6,] 0.007054308 0.003324987 [7,] 0.006932428 0.003266774 [8,] 0.006838894 0.003212682 [9,] 0.006727494 0.003157938 [10,] 0.006620084 0.003105573 [11,] 0.006486799 0.003062020 [12,] 0.006360725 0.003017919 [13,] 0.006258151 0.002975044 [14,] 0.006142128 0.002932265 [15,] 0.006037006 0.002893804 [16,] 0.005933571 0.002854530 [17,] 0.005846824 0.002819702 [18,] 0.005757310 0.002780147 [19,] 0.005674641 0.002738215 [20,] 0.005589567 0.002698052 [21,] 0.005505755 0.002655832 [22,] 0.005447538 0.002616035 [23,] 0.005402345 0.002574503 [24,] 0.005349008 0.002533940 [25,] 0.005272838 0.002493320 [26,] 0.005201531 0.002453994 [27,] 0.005109531 0.002415137 [28,] 0.005044973 0.002380408 [29,] 0.005030640 0.002352465 [30,] 0.004998547 0.002322870 [31,] 0.004947020 0.002293513 [32,] 0.004882920 0.002260166 [33,] 0.004793133 0.002224564 [34,] 0.004727957 0.002192368 [35,] 0.004643047 0.002157118 [36,] 0.004558362 0.002115368 [37,] 0.004495405 0.002072734 [38,] 0.004437511 0.002021947 [39,] 0.004350175 0.001975713 [40,] 0.004223406 0.001933234 [41,] 0.004143031 0.001896464 [42,] 0.004059630 0.001850596 [43,] 0.003958364 0.001798402 > (midr <- midrange(x)) [1] 0.009247167 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 0.004356542 0.004882920 0.004882920 0.004882920 0.004882920 0.004436211 [7] 0.004882920 0.004882920 > postscript(file="/var/www/html/rcomp/tmp/12f3u1197041670.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/2pphw1197041670.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/3r0zc1197041671.tab") > > system("convert tmp/12f3u1197041670.ps tmp/12f3u1197041670.png") > system("convert tmp/2pphw1197041670.ps tmp/2pphw1197041670.png") > > > proc.time() user system elapsed 1.216 0.355 1.375