R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-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(237.588
+ ,164.083
+ ,278.261
+ ,220.36
+ ,253.967
+ ,422.31
+ ,136.921
+ ,143.495
+ ,189.785
+ ,219.529
+ ,217.761
+ ,221.754
+ ,159.854
+ ,209.464
+ ,174.283
+ ,154.55
+ ,153.024
+ ,162.49
+ ,154.462
+ ,249.671
+ ,259.473
+ ,155.337
+ ,151.289
+ ,276.614
+ ,188.214
+ ,181.098
+ ,240.898
+ ,244.551
+ ,250.238
+ ,183.129
+ ,310.331
+ ,281.942
+ ,230.343
+ ,161.563
+ ,392.527
+ ,1077.414
+ ,248.275
+ ,557.386
+ ,731.874
+ ,301.429
+ ,226.36
+ ,215.018
+ ,157.672
+ ,219.118
+ ,213.019
+ ,390.642
+ ,157.124
+ ,227.652
+ ,239.266
+ ,506.343
+ ,149.219
+ ,213.351
+ ,174.517
+ ,172.531
+ ,320.656
+ ,305.011
+ ,266.495
+ ,361.511
+ ,361.019
+ ,382.187
+ ,196.763
+ ,273.212
+ ,186.397
+ ,294.205
+ ,364.685
+ ,230.501
+ ,217.51
+ ,262.297
+ ,169.246
+ ,260.428
+ ,348.187
+ ,512.937
+ ,164.496
+ ,111.187
+ ,169.999
+ ,240.187
+ ,187.158
+ ,194.096
+ ,265.846
+ ,283.319
+ ,356.938
+ ,240.802
+ ,326.662
+ ,249.266
+ ,277.368
+ ,394.618
+ ,235.686
+ ,227.641
+ ,159.593
+ ,268.866
+ ,206.466
+ ,233.064
+ ,133.824
+ ,486.783
+ ,228.859
+ ,155.238
+ ,2042.451
+ ,205.218
+ ,373.648
+ ,229.151
+ ,199.156
+ ,234.41
+ ,56.519
+ ,289.239
+ ,199.227
+ ,274.513
+ ,174.499
+ ,217.714
+ ,239.717
+ ,241.529
+ ,155.561
+ ,204.107
+ ,745.97
+ ,241.772
+ ,110.267
+ ,186.58
+ ,227.906
+ ,197.518
+ ,254.094
+ ,173.942
+ ,294.42
+ ,211.924
+ ,262.479
+ ,193.495
+ ,165.972
+ ,237.352
+ ,205.814
+ ,227.526
+ ,250.439
+ ,470.849
+ ,176.469
+ ,298.691
+ ,193.922
+ ,212.422
+ ,203.284
+ ,240.56
+ ,445.327
+ ,248.984
+ ,174.44
+ ,165.024
+ ,249.681
+ ,238.312
+ ,250.437
+ ,174.75
+ ,4941.633
+ ,138.936
+ ,203.181
+ ,187.747
+ ,270.95
+ ,307.688
+ ,184.477
+ ,230.916
+ ,187.286
+ ,169.376
+ ,182.838
+ ,176.081
+ ,248.056
+ ,235.24
+ ,76.347)
> 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] 288.8592
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 32.91641
> (armose <- arm / armse)
[1] 8.775537
> (geo <- geomean(x))
[1] 236.9969
> (har <- harmean(x))
[1] 216.6176
> (qua <- quamean(x))
[1] 504.6098
> (win <- winmean(x))
          [,1]      [,2]
 [1,] 270.7500 18.486283
 [2,] 259.0379 11.909047
 [3,] 252.8016  9.480633
 [4,] 253.0164  9.314544
 [5,] 247.6268  7.696223
 [6,] 246.0255  7.280462
 [7,] 245.9359  7.195143
 [8,] 245.2398  6.948509
 [9,] 244.4550  6.743168
[10,] 242.9590  6.402185
[11,] 241.4661  6.083949
[12,] 239.3828  5.706285
[13,] 239.2681  5.671664
[14,] 239.1108  5.642547
[15,] 238.3343  5.506801
[16,] 237.6323  5.351661
[17,] 236.7326  5.194234
[18,] 236.5907  5.117636
[19,] 236.5631  5.105591
[20,] 236.2648  5.006460
[21,] 235.2314  4.820930
[22,] 232.4735  4.374907
[23,] 231.6645  4.252212
[24,] 230.1857  4.043045
[25,] 229.9192  3.976062
[26,] 230.0168  3.868234
[27,] 229.4306  3.792158
[28,] 229.0582  3.723647
[29,] 228.7410  3.586293
[30,] 228.9666  3.555342
[31,] 228.0649  3.438027
[32,] 226.9051  3.303343
[33,] 226.6315  3.271372
[34,] 225.8482  3.187898
[35,] 225.7029  3.162279
[36,] 225.8336  3.114039
[37,] 225.4350  3.054803
[38,] 226.2303  2.912717
[39,] 226.1023  2.814724
[40,] 225.6512  2.755125
[41,] 225.3874  2.660991
[42,] 225.7232  2.595223
[43,] 224.8621  2.502569
[44,] 224.9717  2.482323
[45,] 224.4789  2.428872
[46,] 224.3360  2.390153
[47,] 222.8840  2.232240
[48,] 223.3200  2.182466
[49,] 223.3761  1.976649
[50,] 223.5097  1.963768
[51,] 223.5017  1.952913
[52,] 224.1917  1.854698
[53,] 224.4401  1.830919
> (tri <- trimean(x))
          [,1]      [,2]
 [1,] 260.7036 14.838529
 [2,] 250.3978  9.529401
 [3,] 245.9084  7.925918
 [4,] 243.4889  7.253988
 [5,] 240.9471  6.530085
 [6,] 239.5022  6.215314
 [7,] 239.5022  5.969244
 [8,] 237.0987  5.713851
 [9,] 235.9511  5.479544
[10,] 234.8703  5.257395
[11,] 233.9315  5.070789
[12,] 233.1248  4.916061
[13,] 232.5013  4.802355
[14,] 232.5013  4.681626
[15,] 231.2320  4.551990
[16,] 230.6392  4.426775
[17,] 230.0833  4.308177
[18,] 229.5777  4.196728
[19,] 229.0657  4.082125
[20,] 228.5385  3.955695
[21,] 228.0135  3.826457
[22,] 227.5383  3.704957
[23,] 227.2226  3.622668
[24,] 226.9460  3.544740
[25,] 226.7491  3.480801
[26,] 226.5606  3.416388
[27,] 226.3593  3.355287
[28,] 226.3593  3.294815
[29,] 226.1837  3.234015
[30,] 225.8716  3.179710
[31,] 225.7025  3.121325
[32,] 225.5749  3.067755
[33,] 225.5038  3.021113
[34,] 225.4441  2.971263
[35,] 225.4229  2.923116
[36,] 225.4083  2.870625
[37,] 225.3862  2.815525
[38,] 225.3837  2.758625
[39,] 225.3399  2.708720
[40,] 225.3006  2.661773
[41,] 225.2825  2.613833
[42,] 225.2770  2.568626
[43,] 225.2539  2.523123
[44,] 225.2743  2.480470
[45,] 225.2902  2.432180
[46,] 225.3329  2.381351
[47,] 225.3860  2.325547
[48,] 225.5203  2.279214
[49,] 225.6398  2.230186
[50,] 225.7643  2.199480
[51,] 225.8901  2.162788
[52,] 226.0255  2.118253
[53,] 226.1312  2.078809
> (midr <- midrange(x))
[1] 2499.076
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 224.7698 225.3399 225.3399 225.3399 225.3006 224.7698 225.3399 225.3399
> postscript(file="/var/www/rcomp/tmp/1p3tc1291978947.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/www/rcomp/tmp/2p3tc1291978947.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/3ro2t1291978947.tab") 
> 
> try(system("convert tmp/1p3tc1291978947.ps tmp/1p3tc1291978947.png",intern=TRUE))
character(0)
> try(system("convert tmp/2p3tc1291978947.ps tmp/2p3tc1291978947.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.940   0.430   1.368