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.

  Natural language support but running in an English locale

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(426.113
+ ,383.703
+ ,232.444
+ ,70.939
+ ,226.731
+ ,947.293
+ ,611.281
+ ,158.047
+ ,33.999
+ ,37.028
+ ,388.3
+ ,506.652
+ ,392.25
+ ,180.818
+ ,198.296
+ ,217.465
+ ,275.562
+ ,1030.944
+ ,57.47
+ ,136.452
+ ,556.277
+ ,213.361
+ ,274.482
+ ,220.553
+ ,236.71
+ ,260.642
+ ,2763.544
+ ,213.923
+ ,169.861
+ ,403.064
+ ,449.594
+ ,406.167
+ ,206.893
+ ,156.187
+ ,257.102
+ ,62.156
+ ,662.883
+ ,251.422
+ ,171.328
+ ,350.089
+ ,221.588
+ ,4.813
+ ,183.186
+ ,190.379
+ ,223.166
+ ,232.669
+ ,356.725
+ ,109.215
+ ,475.834
+ ,315.955
+ ,694.87
+ ,8.95
+ ,278.741
+ ,308.16
+ ,207.533
+ ,192.797
+ ,601.162
+ ,289.714
+ ,293.671
+ ,386.688
+ ,699.645
+ ,85.094
+ ,131.812
+ ,645.285
+ ,197.549
+ ,308.174
+ ,86.58
+ ,242.205
+ ,238.502
+ ,187.881
+ ,140.321
+ ,440.31
+ ,421.403
+ ,218.761
+ ,1305.923
+ ,137.55
+ ,262.517
+ ,348.821
+ ,150.034
+ ,64.016
+ ,261.596
+ ,259.7
+ ,171.26
+ ,203.077
+ ,249.148
+ ,211.655
+ ,252.64
+ ,438.555
+ ,239.89
+ ,401.915
+ ,216.886
+ ,184.641
+ ,380.155
+ ,653.641
+ ,313.906
+ ,366.936
+ ,236.302
+ ,229.641
+ ,235.577
+ ,103.898
+ ,263.906
+ ,241.171
+ ,216.548
+ ,295.281
+ ,193.299
+ ,204.386
+ ,257.567
+ ,136.813
+ ,240.755
+ ,59.609
+ ,213.511
+ ,380.531
+ ,242.344
+ ,250.407
+ ,183.613
+ ,191.835
+ ,266.793
+ ,246.542
+ ,330.563
+ ,403.556
+ ,208.108
+ ,324.04
+ ,308.532
+ ,199.297
+ ,200.156
+ ,262.875
+ ,287.069
+ ,190.157
+ ,199.746
+ ,265.777
+ ,435.956
+ ,72.844
+ ,756.46
+ ,206.771
+ ,4202.361
+ ,401.422
+ ,216.046
+ ,39.047
+ ,441.437)
> 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] 333.3408
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 36.88561
> (armose <- arm / armse)
[1] 9.037152
> (geo <- geomean(x))
[1] 240.5751
> (har <- harmean(x))
[1] 137.5824
> (qua <- quamean(x))
[1] 546.6916
> (win <- winmean(x))
          [,1]      [,2]
 [1,] 323.0194 29.784872
 [2,] 302.4068 19.031524
 [3,] 296.5374 16.947215
 [4,] 294.1883 16.210742
 [5,] 287.9865 14.306138
 [6,] 285.6264 13.735977
 [7,] 285.5142 13.668381
 [8,] 283.7802 13.265087
 [9,] 283.6301 13.088525
[10,] 283.1660 12.949907
[11,] 281.4445 12.302453
[12,] 280.6992 12.119025
[13,] 278.1209 11.150695
[14,] 273.6583 10.223038
[15,] 272.7711  9.413039
[16,] 270.2848  8.897906
[17,] 269.3313  8.749120
[18,] 269.2808  8.717882
[19,] 269.4197  8.642583
[20,] 270.4433  8.441598
[21,] 269.8858  8.137356
[22,] 269.4347  8.004470
[23,] 268.8685  7.470388
[24,] 268.6592  7.387401
[25,] 268.5830  7.374532
[26,] 270.1432  7.181160
[27,] 270.5074  7.127270
[28,] 268.7458  6.876354
[29,] 268.1362  6.750336
[30,] 268.4875  6.643384
[31,] 268.3294  6.513352
[32,] 267.6503  6.415477
[33,] 267.9067  6.374058
[34,] 264.9086  5.946035
[35,] 262.4638  5.618614
[36,] 261.8459  5.313319
[37,] 261.7072  5.255326
[38,] 256.9895  4.643955
[39,] 255.2852  4.424404
[40,] 253.0766  4.154919
[41,] 253.3338  4.009802
[42,] 252.1055  3.797501
[43,] 252.7326  3.718922
[44,] 252.7668  3.714986
[45,] 248.8045  3.254603
[46,] 248.4620  3.182122
> (tri <- trimean(x))
          [,1]      [,2]
 [1,] 307.4978 24.219556
 [2,] 291.5163 16.244255
 [3,] 285.8254 14.466021
 [4,] 282.0367 13.393383
 [5,] 278.7633 12.431731
 [6,] 276.7444 11.921432
 [7,] 275.0983 11.494407
 [8,] 273.4168 11.031020
 [9,] 271.9286 10.595819
[10,] 270.4099 10.140042
[11,] 268.8945  9.650435
[12,] 267.5155  9.212893
[13,] 266.1640  8.747759
[14,] 265.0123  8.381934
[15,] 264.2247  8.114379
[16,] 263.4846  7.927191
[17,] 262.9219  7.785635
[18,] 262.4131  7.645600
[19,] 261.8880  7.491674
[20,] 261.3315  7.326949
[21,] 260.6786  7.164022
[22,] 260.0371  7.016183
[23,] 259.3987  6.864044
[24,] 258.7698  6.753309
[25,] 258.1262  6.635722
[26,] 257.4579  6.501755
[27,] 256.6601  6.367477
[28,] 255.8012  6.217172
[29,] 255.0079  6.075572
[30,] 254.2113  5.926427
[31,] 253.3523  5.763902
[32,] 252.4569  5.589137
[33,] 251.5528  5.396337
[34,] 250.5826  5.170835
[35,] 249.7338  4.971327
[36,] 248.9793  4.785975
[37,] 248.2150  4.609307
[38,] 247.4104  4.403133
[39,] 246.8360  4.261869
[40,] 246.3256  4.127045
[41,] 245.9140  4.008750
[42,] 245.4566  3.884973
[43,] 245.0415  3.769089
[44,] 244.5540  3.634764
[45,] 244.0245  3.463337
[46,] 243.7103  3.348964
> (midr <- midrange(x))
[1] 2103.587
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 248.9204 250.5826 250.5826 250.5826 249.7338 248.9204 250.5826 250.5826
> postscript(file="/var/www/html/freestat/rcomp/tmp/13ugb1286218919.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/2wlxv1286218919.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/3z76d1286218919.tab") 
> try(system("convert tmp/13ugb1286218919.ps tmp/13ugb1286218919.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wlxv1286218919.ps tmp/2wlxv1286218919.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.219   0.543   1.339