R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: i686-pc-linux-gnu (32-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(0.0559998852152633
+ ,-0.493896873364634
+ ,-0.493896873364634
+ ,-0.987774567512833
+ ,9.93823000176787
+ ,-0.586194140197587
+ ,-7.68751056589317
+ ,-3.33966294852904
+ ,0.429813645031793
+ ,-1.22306236050721
+ ,-0.387148988105899
+ ,1.0240939117167
+ ,-3.1790261312384
+ ,-5.0454536432582
+ ,-4.0137177585503
+ ,1.44819243280776
+ ,0.841409399737355
+ ,2.17969801643555
+ ,-1.73266657002515
+ ,-8.06571228807107
+ ,8.60981291854194
+ ,-0.334926871345123
+ ,4.75053530730328
+ ,3.61827550903922
+ ,4.19213041938715
+ ,-3.79054822946422
+ ,9.70309656335496
+ ,-1.30481000432226
+ ,-0.239045671970004
+ ,-4.18403094181251
+ ,-3.35133824921843
+ ,2.80361259428111
+ ,2.26146839403967
+ ,2.49924097981172
+ ,-2.60790853259369
+ ,1.07120014090196
+ ,-6.63714848357267
+ ,1.14652131023368
+ ,-1.23746151713059
+ ,1.54755663980393
+ ,11.5631756773691
+ ,3.24139827235349
+ ,-0.701479351284132
+ ,10.2563339658872
+ ,-2.94489100586251
+ ,-5.78748839028522
+ ,6.79431478166191
+ ,5.9550320151792
+ ,7.99483304039295
+ ,-7.04636839223017
+ ,-5.51175848427593
+ ,-3.38566124618367
+ ,3.18364393412584
+ ,0.906160332366917
+ ,7.17349376234177
+ ,7.47708219897042
+ ,-5.86046612170084
+ ,-4.08054017321996
+ ,3.19332756653989
+ ,0.319595007657352
+ ,1.37314374171947
+ ,0.459430792130583
+ ,-4.31529513951182
+ ,9.8377414934571
+ ,2.65447899718868
+ ,5.51664743247295
+ ,3.89578300556218
+ ,0.110260870088769
+ ,0.905004879136158
+ ,-2.69130528473603
+ ,9.15761451252449
+ ,-4.67587915706826
+ ,2.02828775800485
+ ,-6.08812584767335
+ ,-0.848867752979329
+ ,3.22884771065786
+ ,-4.76216620747938
+ ,10.0643639249429
+ ,-1.80365481397992
+ ,-4.041705411162
+ ,4.53647207420368
+ ,1.27722948511292
+ ,2.71494904872247
+ ,-2.98625229115164
+ ,-7.91791855375443
+ ,-0.87819280547091
+ ,8.48347072627108
+ ,-6.28135915023324
+ ,-3.67782642483031
+ ,3.79098609354026
+ ,-9.35494494875011
+ ,5.42924781064487
+ ,-6.35112635227124
+ ,1.64151446618105
+ ,8.41860163489749
+ ,-6.11616812572276
+ ,-3.89690540156151
+ ,5.94210917809354
+ ,6.57477404609063
+ ,2.53897589056262
+ ,-4.46837465429453
+ ,-5.48322577209247
+ ,2.92641420148313
+ ,-1.06268550612349
+ ,-2.78613020257466
+ ,1.99145230799011
+ ,5.06497199204414
+ ,6.04930091019776
+ ,-1.6778355872977
+ ,-4.50123046727257
+ ,4.39563205247514
+ ,0.351093295039824
+ ,-0.078743595656535
+ ,2.63210831940019
+ ,-1.96551715036875
+ ,2.80719539904876
+ ,2.53659023294771
+ ,-10.4791047436627
+ ,1.80942984068822
+ ,-1.36718444180429)
> 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] 0.421034
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 0.4378202
> (armose <- arm / armse)
[1] 0.9616596
> (geo <- geomean(x))
[1] NaN
Warning message:
In log(x) : NaNs produced
> (har <- harmean(x))
[1] 10.85612
> (qua <- quamean(x))
[1] 4.794577
> (win <- winmean(x))
            [,1]      [,2]
 [1,] 0.41951166 0.4337500
 [2,] 0.43779937 0.4293197
 [3,] 0.43834087 0.4281181
 [4,] 0.44267152 0.4262502
 [5,] 0.46377557 0.4210696
 [6,] 0.45696246 0.4131305
 [7,] 0.44169199 0.4052160
 [8,] 0.43792033 0.4031430
 [9,] 0.44544447 0.4006078
[10,] 0.41246728 0.3944793
[11,] 0.38587560 0.3840915
[12,] 0.36281452 0.3784435
[13,] 0.35160754 0.3682783
[14,] 0.32932327 0.3641035
[15,] 0.31836065 0.3475937
[16,] 0.34356312 0.3410282
[17,] 0.35395638 0.3392444
[18,] 0.31633442 0.3273119
[19,] 0.30769832 0.3248283
[20,] 0.27249894 0.3137639
[21,] 0.24044375 0.3039844
[22,] 0.22017213 0.2968241
[23,] 0.20062112 0.2926422
[24,] 0.16551833 0.2871086
[25,] 0.12811519 0.2769653
[26,] 0.12845325 0.2715763
[27,] 0.11495577 0.2641402
[28,] 0.09518963 0.2463521
[29,] 0.10045130 0.2450432
[30,] 0.09449009 0.2437470
[31,] 0.13348633 0.2386014
[32,] 0.11629809 0.2253473
[33,] 0.09488727 0.2206192
[34,] 0.13885437 0.2153422
[35,] 0.14065144 0.2095132
[36,] 0.14752945 0.2048083
[37,] 0.33870250 0.1823365
[38,] 0.36046697 0.1738208
[39,] 0.38276281 0.1713352
[40,] 0.38859005 0.1681289
> (tri <- trimean(x))
           [,1]      [,2]
 [1,] 0.4189831 0.4250523
 [2,] 0.4184364 0.4154097
 [3,] 0.4082454 0.4072948
 [4,] 0.3974970 0.3987618
 [5,] 0.3851766 0.3898549
 [6,] 0.3677102 0.3813130
 [7,] 0.3508701 0.3736464
 [8,] 0.3358995 0.3667377
 [9,] 0.3208964 0.3593785
[10,] 0.3042900 0.3515474
[11,] 0.2910438 0.3437813
[12,] 0.2802675 0.3367122
[13,] 0.2714859 0.3295900
[14,] 0.2634469 0.3230366
[15,] 0.2571730 0.3161843
[16,] 0.2516105 0.3107089
[17,] 0.2435913 0.3053194
[18,] 0.2343170 0.2993482
[19,] 0.2276489 0.2940810
[20,] 0.2213292 0.2883036
[21,] 0.2173931 0.2830748
[22,] 0.2156599 0.2782611
[23,] 0.2153273 0.2735467
[24,] 0.2163930 0.2685316
[25,] 0.2200269 0.2633284
[26,] 0.2265148 0.2585067
[27,] 0.2333722 0.2534667
[28,] 0.2415956 0.2484342
[29,] 0.2517158 0.2448254
[30,] 0.2621479 0.2405498
[31,] 0.2737105 0.2354455
[32,] 0.2834034 0.2299579
[33,] 0.2950079 0.2252329
[34,] 0.3090024 0.2200119
[35,] 0.3210128 0.2143433
[36,] 0.3338958 0.2081533
[37,] 0.3474006 0.2010838
[38,] 0.3480417 0.1966021
[39,] 0.3471075 0.1923332
[40,] 0.3443648 0.1870498
> (midr <- midrange(x))
[1] 0.5420355
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 0.2029104 0.2621479 0.2029104 0.2621479 0.2621479 0.2029104 0.2621479
[8] 0.2517158
> postscript(file="/var/fisher/rcomp/tmp/1b3pd1387556643.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/fisher/rcomp/tmp/24wwo1387556643.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/3z0h41387556644.tab") 
> 
> try(system("convert tmp/1b3pd1387556643.ps tmp/1b3pd1387556643.png",intern=TRUE))
character(0)
> try(system("convert tmp/24wwo1387556643.ps tmp/24wwo1387556643.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  3.450   0.814   4.231