R version 2.8.0 (2008-10-20) Copyright (C) 2008 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(10.0699948358949 + ,66.1571577530624 + ,-151.075300538186 + ,-248.831522488301 + ,-639.593888718853 + ,76.1920145357235 + ,192.845671000015 + ,185.456059757973 + ,545.777877412695 + ,204.807985464276 + ,-188.598125817465 + ,-361.597229103207 + ,237.791037067447 + ,88.6706048611213 + ,-1.67573867458668 + ,-101.324261325413 + ,266.354817734124 + ,664.11529839286 + ,-127.829170960315 + ,140.359503190461 + ,-46.8833565357527 + ,-148.337869337294 + ,-690.093664539694 + ,10.9872887023776 + ,297.880016150804 + ,15.5158275040067 + ,-127.162354841272 + ,224.150091900778 + ,315.072479043656 + ,-54.3933999841192 + ,-609.217437130968 + ,-497.593440361128 + ,187.700712912703 + ,20.6527597500335 + ,167.152983928598 + ,419.171277396816 + ,485.910124202384 + ,-448.179303238101 + ,225.321369297593 + ,-122.105927480186 + ,-23.5350176864831 + ,-319.802851650811 + ,345.431533876987 + ,251.546832269847 + ,-532.83385641665 + ,33.5562031825193 + ,10.5260951309392 + ,234.075819428605 + ,65.2437564404718 + ,144.410348380949 + ,44.3891628849124 + ,-259.491301623020 + ,-242.033896793660 + ,45.4061112817417 + ,132.055082289697 + ,124.469667769854 + ,47.1454064444406 + ,636.29504998809 + ,-21.2305967857192 + ,96.1454064444406 + ,410.065349916627 + ,611.106824194445 + ,440.397637083328 + ,699.885598321398 + ,635.548177341234 + ,314.740191876957 + ,-168.489956551633 + ,1003.25781280948 + ,997.808433821407 + ,-1888.88470160714 + ,705.066246630886 + ,-378.446688845301 + ,-402.208761047506 + ,-322.459848500052 + ,727.362843575409 + ,-232.001345071387 + ,302.935546797629 + ,659.266287007895 + ,-365.765614472202 + ,199.675738674587 + ,-355.592543646871 + ,-157.590548333308 + ,-157.422611321443 + ,125.999551642872 + ,790.418374222234 + ,612.816011305564 + ,649.927520956342 + ,-70.3788952579816 + ,1555.07806121425 + ,649.009819269802 + ,-47.7477693611145 + ,-1784.73615666280 + ,-262.220777515919 + ,-43.6752903174565 + ,1027.97143890485 + ,-1076.38916288491 + ,411.905887103178 + ,1036.56133699599 + ,-275.228601472158) > 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] 79.5248 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 50.96049 > (armose <- arm / armse) [1] 1.560519 > (geo <- geomean(x)) [1] NaN Warning message: In log(x) : NaNs produced > (har <- harmean(x)) [1] -480.6038 > (qua <- quamean(x)) [1] 510.7127 > (win <- winmean(x)) [,1] [,2] [1,] 75.33926 49.25081 [2,] 89.47577 44.50200 [3,] 100.43280 41.62535 [4,] 102.25302 41.19132 [5,] 93.31293 38.78002 [6,] 94.12070 37.27211 [7,] 95.03591 36.58142 [8,] 98.61034 35.86920 [9,] 99.53763 34.68718 [10,] 101.44804 34.25709 [11,] 101.81941 33.88534 [12,] 102.21343 33.79603 [13,] 101.33231 33.41158 [14,] 105.91212 32.75453 [15,] 102.87043 32.13762 [16,] 109.79812 31.15450 [17,] 100.81368 29.07288 [18,] 90.42491 27.35456 [19,] 83.73600 25.82607 [20,] 80.82111 25.04775 [21,] 81.40808 24.55718 [22,] 90.64423 23.28119 [23,] 80.29989 20.67064 [24,] 75.58240 19.40061 [25,] 75.54090 19.38484 [26,] 74.10766 18.79290 [27,] 73.47545 18.53082 [28,] 70.35967 16.74167 [29,] 66.21731 16.20570 [30,] 63.58114 15.53721 [31,] 68.92518 14.61990 [32,] 76.09800 13.12771 [33,] 81.03607 12.48376 > (tri <- trimean(x)) [,1] [,2] [1,] 84.60579 45.33834 [2,] 94.26249 40.61662 [3,] 96.81026 38.20582 [4,] 95.49660 36.72829 [5,] 93.61771 35.18012 [6,] 93.68707 34.12655 [7,] 93.60289 33.30451 [8,] 93.35871 32.51667 [9,] 92.55638 31.75676 [10,] 91.58431 31.10921 [11,] 90.31611 30.43365 [12,] 88.93572 29.70984 [13,] 87.43515 28.87416 [14,] 85.94456 27.95634 [15,] 83.89819 26.97917 [16,] 82.02929 25.92577 [17,] 79.38591 24.82738 [18,] 77.40519 23.89254 [19,] 76.23128 23.07800 [20,] 75.56851 22.36059 [21,] 75.11236 21.62016 [22,] 74.57273 20.79461 [23,] 73.20817 19.99966 [24,] 72.60963 19.50184 [25,] 72.35937 19.10005 [26,] 72.09131 18.58159 [27,] 72.09131 18.02624 [28,] 71.78812 17.35276 [29,] 71.91131 16.85874 [30,] 72.40972 16.30940 [31,] 73.19714 15.71606 [32,] 73.58693 15.13668 [33,] 73.35152 14.72283 > (midr <- midrange(x)) [1] -166.9033 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 67.76037 72.60963 72.60963 72.60963 72.35937 67.76037 72.60963 72.60963 > postscript(file="/var/www/html/freestat/rcomp/tmp/1gp7c1229288880.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/freestat/rcomp/tmp/2957b1229288880.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 > > #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/3cv151229288880.tab") > > system("convert tmp/1gp7c1229288880.ps tmp/1gp7c1229288880.png") > system("convert tmp/2957b1229288880.ps tmp/2957b1229288880.png") > > > proc.time() user system elapsed 1.006 0.502 1.110