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(145.3,143.6,142.8,155.9,156.2,149.8,152.7,155.5,159.3,143,141.4,142.8,146.4,152.3,164.3,168,171.3,162.7,150.2,142.5,138.2,138,145.1,138.4,131.8,130.8,126.3,123,124,120.8,122.1,106.5,104.3,108.7,113.8,112.5,106.1,98.4,96,99.3,97.5,95.3,88,94.7,99.4,98.9,96.4,95.3,99.5,101.6,103.9,106.6,108.3,102,93.8,91.6,97.7,94.8,98,103.8,97.8,91.2,89.3,87.5,90.4,94.2,102.2,101.3,96,90.8,93.2,90.9,91.1,90.2,94.3,96,99,103.3,113.1,112.8,112.1,107.4,111,110.5,110.8,112.4,111.5,116.2,122.5,121.3,113.9,110.7,120.8,141.1,147.4,148,158.1,165,187,190.3,182.4,168.8,151.2,120.1,112.5,106.2,107.1,108.5,106.5,108.3,125.6,124,127.2,136.9,135.8,124.3,115.4,113.6,114.4,118.4,117,116.5,115.4,113.6,117.4,116.9,116.4,111.1,110.2,118.9,131.8,130.6,138.3,148.4,148.7,144.3,152.5,162.9,167.2,166.5,185.6)
> 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] 122.5333
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 2.101339
> (armose <- arm / armse)
[1] 58.31204
> (geo <- geomean(x))
[1] 120.1581
> (har <- harmean(x))
[1] 117.9416
> (qua <- quamean(x))
[1] 125.0304
> (win <- winmean(x))
          [,1]      [,2]
 [1,] 122.5135 2.0956510
 [2,] 122.5121 2.0891953
 [3,] 122.4631 2.0726678
 [4,] 122.1539 2.0119651
 [5,] 122.0794 1.9952221
 [6,] 122.0496 1.9890902
 [7,] 122.0199 1.9814705
 [8,] 121.9858 1.9744060
 [9,] 121.9156 1.9563303
[10,] 121.9794 1.9362448
[11,] 121.9170 1.9144274
[12,] 121.9340 1.9082747
[13,] 121.6298 1.8607147
[14,] 121.5504 1.8394899
[15,] 121.3589 1.8101361
[16,] 121.3816 1.7995538
[17,] 121.3333 1.7929677
[18,] 121.0652 1.7361202
[19,] 121.0383 1.7326186
[20,] 121.0099 1.7289452
[21,] 120.9057 1.7017974
[22,] 120.9213 1.6646217
[23,] 120.8887 1.6531739
[24,] 120.7184 1.6282783
[25,] 120.7007 1.6182027
[26,] 120.7007 1.6018412
[27,] 120.6816 1.5784223
[28,] 120.5028 1.5526253
[29,] 120.3383 1.5197723
[30,] 120.3170 1.5126836
[31,] 120.1631 1.4900410
[32,] 120.4128 1.4319420
[33,] 120.3426 1.4090737
[34,] 120.3908 1.3944106
[35,] 120.4404 1.3897450
[36,] 120.6447 1.3549221
[37,] 120.4872 1.3098628
[38,] 120.4333 1.2982038
[39,] 119.7972 1.2043849
[40,] 120.2794 1.1553696
[41,] 120.2794 1.1495845
[42,] 120.3092 1.1351606
[43,] 119.9738 1.0980173
[44,] 119.6617 1.0576659
[45,] 118.5447 0.9070267
[46,] 118.6426 0.8982411
[47,] 118.6092 0.8368404
> (tri <- trimean(x))
          [,1]      [,2]
 [1,] 122.2978 2.0591940
 [2,] 122.0759 2.0192237
 [3,] 121.8481 1.9790363
 [4,] 121.6308 1.9415092
 [5,] 121.4901 1.9192693
 [6,] 121.3612 1.8990602
 [7,] 121.2339 1.8781172
 [8,] 121.1072 1.8564848
 [9,] 120.9813 1.8338344
[10,] 120.8603 1.8117639
[11,] 120.7277 1.7903096
[12,] 120.5974 1.7696301
[13,] 120.4609 1.7474042
[14,] 120.3487 1.7290388
[15,] 120.2396 1.7111602
[16,] 120.1431 1.6947278
[17,] 120.0411 1.6774233
[18,] 119.9390 1.6586710
[19,] 119.8534 1.6440270
[20,] 119.7663 1.6277399
[21,] 119.6778 1.6096068
[22,] 119.5928 1.5919709
[23,] 119.5032 1.5757971
[24,] 119.4118 1.5584437
[25,] 119.3275 1.5412458
[26,] 119.2404 1.5224778
[27,] 119.1494 1.5025567
[28,] 119.0553 1.4820613
[29,] 118.9675 1.4612470
[30,] 118.8852 1.4407914
[31,] 118.8000 1.4176657
[32,] 118.7195 1.3933369
[33,] 118.6200 1.3713793
[34,] 118.5192 1.3481605
[35,] 118.4099 1.3222811
[36,] 118.2913 1.2919083
[37,] 118.1537 1.2598779
[38,] 118.0169 1.2278197
[39,] 117.8746 1.1907197
[40,] 117.7607 1.1607923
[41,] 117.6102 1.1303752
[42,] 117.4491 1.0935048
[43,] 117.2745 1.0497626
[44,] 117.1075 1.0025760
[45,] 116.9471 0.9515904
[46,] 116.8449 0.9200620
[47,] 116.7277 0.8808514
> (midr <- midrange(x))
[1] 138.9
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 118.0614 118.7486 118.7486 118.7486 118.7486 118.5192 118.7486 118.7486
> postscript(file="/var/www/html/rcomp/tmp/1o5ya1196431244.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/2doos1196431244.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/356ks1196431245.tab") 
> 
> system("convert tmp/1o5ya1196431244.ps tmp/1o5ya1196431244.png")
> system("convert tmp/2doos1196431244.ps tmp/2doos1196431244.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.252   0.344   1.364