R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(294.371
+ ,155.359
+ ,167.804
+ ,82.893
+ ,103.131
+ ,130.093
+ ,166.796
+ ,72.174
+ ,96.396
+ ,81.056
+ ,138.981
+ ,197.604
+ ,264.405
+ ,119.403
+ ,189.62
+ ,155.779
+ ,134.168
+ ,126.862
+ ,159.104
+ ,163.426
+ ,113.344
+ ,175.645
+ ,147.312
+ ,164.642
+ ,125.814
+ ,138.558
+ ,163.136
+ ,148.763
+ ,178.322
+ ,107.805
+ ,180.363
+ ,179.006
+ ,111.603
+ ,251.469
+ ,123.646
+ ,108.398
+ ,134.505
+ ,129.605
+ ,211.562
+ ,300.239
+ ,329.183
+ ,115.475
+ ,115.108
+ ,101.922
+ ,391.433
+ ,103.906
+ ,112.704
+ ,619.22
+ ,110.688
+ ,95.539
+ ,150.985
+ ,138.107
+ ,177.705
+ ,757.365
+ ,103.671
+ ,81.101
+ ,130.586
+ ,166.608
+ ,201.734
+ ,175.329
+ ,126.372
+ ,107.822
+ ,146.123
+ ,98.815
+ ,219.641
+ ,156
+ ,123.928
+ ,156.405
+ ,104.828
+ ,188.594
+ ,575.698
+ ,564.424
+ ,149.465
+ ,168.474
+ ,280.909
+ ,172.095
+ ,73.594
+ ,128.317
+ ,202.431
+ ,149.787
+ ,243.719
+ ,188.553
+ ,125.634
+ ,166.548
+ ,284.966
+ ,122.213
+ ,172.936
+ ,152.882
+ ,242.736
+ ,107.814
+ ,155.723
+ ,192.582
+ ,126.668
+ ,127.17
+ ,125.77
+ ,144.493
+ ,123.241
+ ,138.431
+ ,150.875
+ ,375.445
+ ,154.611
+ ,83.401
+ ,115.175
+ ,198.463
+ ,520.108
+ ,235.576
+ ,105.44
+ ,152.272
+ ,310.469
+ ,217.089
+ ,229.36
+ ,555.232
+ ,466.984
+ ,127.703
+ ,101.383
+ ,172.705
+ ,79.114
+ ,213.859
+ ,256.736
+ ,161.995
+ ,112.314
+ ,219.779
+ ,188.723
+ ,141.648
+ ,146.564
+ ,113.163
+ ,85.657
+ ,127.958
+ ,114.909
+ ,198.214
+ ,133.215
+ ,131.576
+ ,287.667
+ ,140.275
+ ,137.435
+ ,190.055
+ ,104.39
+ ,178.032
+ ,79.356
+ ,250.571
+ ,123.115
+ ,247.026
+ ,190.38
+ ,219.114
+ ,131.753
+ ,135.922
+ ,855.826
+ ,145.195
+ ,82.337
+ ,141.758
+ ,78.649
+ ,117.482
+ ,90.986
+ ,101.864
+ ,120.33
+ ,137.019)
> 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] 180.7153
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 9.595725
> (armose <- arm / armse)
[1] 18.83289
> (geo <- geomean(x))
[1] 158.7908
> (har <- harmean(x))
[1] 145.1782
> (qua <- quamean(x))
[1] 216.6335
> (win <- winmean(x))
          [,1]     [,2]
 [1,] 180.0932 9.325536
 [2,] 178.3869 8.674452
 [3,] 177.5589 8.408594
 [4,] 177.2761 8.320578
 [5,] 177.0359 8.228556
 [6,] 175.6867 7.836385
 [7,] 173.3584 7.178433
 [8,] 169.5125 6.213738
 [9,] 168.6194 6.002219
[10,] 165.7986 5.365572
[11,] 164.8547 5.077930
[12,] 164.4180 4.902449
[13,] 164.0005 4.809504
[14,] 163.6159 4.685972
[15,] 163.6031 4.620040
[16,] 163.2364 4.545738
[17,] 161.4442 4.254260
[18,] 160.6988 4.105762
[19,] 160.1231 4.004171
[20,] 160.0381 3.984521
[21,] 159.6260 3.909192
[22,] 159.2214 3.836779
[23,] 159.1667 3.807990
[24,] 158.4290 3.621281
[25,] 157.4343 3.486356
[26,] 155.8388 3.279349
[27,] 155.9146 3.266957
[28,] 156.2310 3.216982
[29,] 156.0247 3.154239
[30,] 155.5403 3.064879
[31,] 155.1613 3.002134
[32,] 153.3824 2.772784
[33,] 153.2733 2.752413
[34,] 152.9015 2.640428
[35,] 152.8902 2.630071
[36,] 152.7649 2.613020
[37,] 151.6450 2.476864
[38,] 151.5975 2.374299
[39,] 151.9965 2.321707
[40,] 152.1226 2.288497
[41,] 152.3818 2.219703
[42,] 152.5899 2.194863
[43,] 152.6133 2.190660
[44,] 150.4175 1.943120
[45,] 150.1074 1.897187
[46,] 150.4088 1.833036
[47,] 150.3624 1.820890
[48,] 150.2753 1.809858
[49,] 149.8036 1.731777
[50,] 149.7971 1.713763
[51,] 149.0782 1.633752
[52,] 149.1039 1.617506
> (tri <- trimean(x))
          [,1]     [,2]
 [1,] 177.0363 8.634867
 [2,] 173.8988 7.831849
 [3,] 171.5650 7.326292
 [4,] 169.4591 6.870948
 [5,] 167.3710 6.382349
 [6,] 165.2769 5.846343
 [7,] 163.3709 5.337113
 [8,] 161.7810 4.920859
 [9,] 160.6886 4.668719
[10,] 159.6778 4.427215
[11,] 158.9652 4.274166
[12,] 158.3324 4.150991
[13,] 157.7239 4.040825
[14,] 157.1354 3.932151
[15,] 156.5623 3.828598
[16,] 155.9718 3.722002
[17,] 155.3912 3.613120
[18,] 154.9284 3.530021
[19,] 154.5045 3.455758
[20,] 154.1069 3.385195
[21,] 153.7010 3.308706
[22,] 153.3081 3.232065
[23,] 152.9269 3.154716
[24,] 152.5350 3.071069
[25,] 152.1736 2.998646
[26,] 151.8579 2.932902
[27,] 151.6238 2.882304
[28,] 151.3759 2.826368
[29,] 151.0998 2.768000
[30,] 150.8239 2.708648
[31,] 150.8239 2.651249
[32,] 150.3114 2.592958
[33,] 150.1451 2.552559
[34,] 149.9770 2.508382
[35,] 149.8210 2.469753
[36,] 149.6582 2.426178
[37,] 149.4940 2.377674
[38,] 149.3806 2.337590
[39,] 149.2639 2.302289
[40,] 149.1201 2.266377
[41,] 148.9619 2.227645
[42,] 148.7812 2.189572
[43,] 148.5791 2.146856
[44,] 148.3638 2.096149
[45,] 148.3638 2.070022
[46,] 148.1531 2.043944
[47,] 148.0297 2.019623
[48,] 147.9007 1.990775
[49,] 147.7676 1.956308
[50,] 147.6519 1.925298
[51,] 147.5279 1.888736
[52,] 147.4367 1.856118
> (midr <- midrange(x))
[1] 464
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 148.8616 149.2639 148.8616 149.2639 149.2639 148.8616 149.2639 149.3806
> postscript(file="/var/www/html/rcomp/tmp/120oz1289577954.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/2v95k1289577954.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> 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/3fuwj1289577954.tab") 
> try(system("convert tmp/120oz1289577954.ps tmp/120oz1289577954.png",intern=TRUE))
character(0)
> try(system("convert tmp/2v95k1289577954.ps tmp/2v95k1289577954.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.851   0.384   1.140