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(0.0326799763806121
+ ,-0.948210140865705
+ ,1.31775695277446
+ ,-5.72384201804306
+ ,1.2568079894424
+ ,3.58561706384319
+ ,-3.09576079606216
+ ,-0.70454278049026
+ ,3.59924783422083
+ ,-2.34471000535835
+ ,-2.989620808753
+ ,1.82848026673821
+ ,0.867406044652756
+ ,-5.31350276508493
+ ,-0.929286157076638
+ ,1.29698608341727
+ ,0.556921007230199
+ ,-0.0233100263777321
+ ,2.55596018827231
+ ,-0.361416691444688
+ ,-1.58934021456799
+ ,-1.10837043086146
+ ,1.47786959811465
+ ,1.03724274024881
+ ,2.09592562691563
+ ,-2.52485282024109
+ ,-3.02854928038647
+ ,5.11841760936073
+ ,2.67335286331825
+ ,-0.0337627205147761
+ ,-3.22198681075089
+ ,-4.0913310871613
+ ,3.44922099409475
+ ,2.87616723248858
+ ,-0.209660432960636
+ ,0.285059576872534
+ ,-3.04811897760346
+ ,2.63115167630307
+ ,-0.534608108493555
+ ,0.93987263670119
+ ,1.47647427083872
+ ,-0.692555149561432
+ ,1.57805388587354
+ ,-0.509641129678686
+ ,3.95104867743037
+ ,-2.24260248313829
+ ,2.10977095770074
+ ,4.59638453198225
+ ,-1.05219825518802
+ ,5.79575946283695
+ ,-6.0187205018209
+ ,-3.39863275328804
+ ,8.46155122144953
+ ,1.59779748036101
+ ,5.87630129482488
+ ,-2.88107929115406
+ ,-2.33019974185080
+ ,5.13716011655738
+ ,2.96459580540549
+ ,5.53664737235919
+ ,-2.24253330315628
+ ,-4.31450690329314
+ ,-0.580875723811879
+ ,4.16454557617349
+ ,7.71449516869876
+ ,-4.47932948405506
+ ,1.15119340334905
+ ,5.5937205174593
+ ,1.57354154675426
+ ,-3.72116627293326
+ ,5.98232888324607
+ ,0.882622434724226
+ ,-10.9608895852359
+ ,-0.191972055215800
+ ,4.2994240157116
+ ,4.05159070596540
+ ,-8.77364595669872
+ ,3.49083917384408
+ ,4.07347502448791
+ ,3.19580183065757
+ ,-3.05938743065843
+ ,0.734364167425412
+ ,5.87813977974351
+ ,-5.50859889253316
+ ,4.37469546783051
+ ,6.01893709617511
+ ,8.75096962085911
+ ,-5.74961974823077
+ ,-0.319494662937842
+ ,5.11906598002663
+ ,9.82275942751156
+ ,4.97559831737166
+ ,10.4814225736623
+ ,4.12290361003019
+ ,-0.305554836701702
+ ,-19.8829430553374
+ ,-3.34762924408284
+ ,-11.1932656693866
+ ,-6.72350884105087
+ ,-1.47847004245831
+ ,5.23368621262139
+ ,-5.29192778347316
+ ,-1.82632816244524
+ ,-1.16317956771913
+ ,2.86415985839632
+ ,6.59723839746982
+ ,-10.0779820362446
+ ,3.33232685697513
+ ,-4.79223897645672
+ ,3.73448262822686
+ ,4.12967929947457)
> 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] 0.3961463
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 0.4441451
> (armose <- arm / armse)
[1] 0.8919298
> (geo <- geomean(x))
[1] NaN
Warning message:
In log(x) : NaNs produced
> (har <- harmean(x))
[1] -1.985558
> (qua <- quamean(x))
[1] 4.675047
> (win <- winmean(x))
           [,1]      [,2]
 [1,] 0.4684977 0.4165146
 [2,] 0.4533732 0.4117124
 [3,] 0.4694134 0.4044559
 [4,] 0.4894956 0.3890955
 [5,] 0.5315172 0.3625563
 [6,] 0.5383544 0.3512581
 [7,] 0.5530161 0.3480939
 [8,] 0.5473648 0.3467311
 [9,] 0.5646679 0.3438743
[10,] 0.5749881 0.3400639
[11,] 0.5571043 0.3369693
[12,] 0.6049546 0.3278204
[13,] 0.6061197 0.3176325
[14,] 0.6147337 0.3130256
[15,] 0.6424474 0.3084407
[16,] 0.6957111 0.3011469
[17,] 0.7232349 0.2917705
[18,] 0.6700116 0.2827469
[19,] 0.6535712 0.2752534
[20,] 0.6627522 0.2707006
[21,] 0.6441161 0.2667428
[22,] 0.6394390 0.2656336
[23,] 0.6420900 0.2649547
[24,] 0.6398197 0.2626226
[25,] 0.6593371 0.2589814
[26,] 0.7192272 0.2459813
[27,] 0.7103675 0.2345242
[28,] 0.6799144 0.2301230
[29,] 0.6992390 0.2270029
[30,] 0.6736421 0.2240561
[31,] 0.7782564 0.2092545
[32,] 0.8128781 0.1977058
[33,] 0.8052510 0.1893983
[34,] 0.8310066 0.1708894
[35,] 0.8204058 0.1659064
[36,] 0.8347296 0.1635492
[37,] 0.8057899 0.1527998
> (tri <- trimean(x))
           [,1]      [,2]
 [1,] 0.4896675 0.4011850
 [2,] 0.5116286 0.3838394
 [3,] 0.5424208 0.3670796
 [4,] 0.5686467 0.3512055
 [5,] 0.5903937 0.3385286
 [6,] 0.6035963 0.3317656
 [7,] 0.6035963 0.3268300
 [8,] 0.6265590 0.3219062
 [9,] 0.6383743 0.3165674
[10,] 0.6483638 0.3110125
[11,] 0.6575152 0.3053326
[12,] 0.6691616 0.2993328
[13,] 0.6761488 0.2939313
[14,] 0.6761488 0.2892709
[15,] 0.6900696 0.2845539
[16,] 0.6945305 0.2797424
[17,] 0.6944241 0.2751818
[18,] 0.6919159 0.2711540
[19,] 0.6937662 0.2676416
[20,] 0.6970736 0.2644666
[21,] 0.6998342 0.2612563
[22,] 0.7042299 0.2578985
[23,] 0.7092591 0.2539585
[24,] 0.7144046 0.2492555
[25,] 0.7200596 0.2438691
[26,] 0.7246292 0.2378549
[27,] 0.7250338 0.2325559
[28,] 0.7250338 0.2278212
[29,] 0.7295869 0.2225934
[30,] 0.7318646 0.2165348
[31,] 0.7318646 0.2093710
[32,] 0.7330616 0.2030927
[33,] 0.7269090 0.1972332
[34,] 0.7207808 0.1912584
[35,] 0.7120039 0.1872066
[36,] 0.7031888 0.1827137
[37,] 0.6922270 0.1769781
> (midr <- midrange(x))
[1] -4.70076
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 0.6712937 0.7250338 0.7250338 0.7250338 0.7261301 0.6712937 0.7250338
[8] 0.7250338
> postscript(file="/var/www/html/rcomp/tmp/1bi1q1262184436.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/2xlry1262184436.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/3z8tb1262184436.tab") 
> 
> try(system("convert tmp/1bi1q1262184436.ps tmp/1bi1q1262184436.png",intern=TRUE))
character(0)
> try(system("convert tmp/2xlry1262184436.ps tmp/2xlry1262184436.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.730   0.316   0.834