R version 2.12.1 (2010-12-16)
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(0.0218,0.0290,0.0263,0.0267,0.0181,0.0133,0.0088,0.0128,0.0126,0.0126,0.0129,0.0110,0.0137,0.0121,0.0174,0.0176,0.0148,0.0104,0.0162,0.0149,0.0179,0.0180,0.0158,0.0186,0.0174,0.0159,0.0126,0.0113,0.0192,0.0261,0.0226,0.0241,0.0226,0.0203,0.0286,0.0255,0.0227,0.0226,0.0257,0.0307,0.0276,0.0251,0.0287,0.0314,0.0311,0.0316,0.0247,0.0257,0.0289,0.0263,0.0238,0.0169,0.0196,0.0219,0.0187,0.0160,0.0163,0.0122,0.0121,0.0149,0.0164,0.0166,0.0177,0.0182,0.0178,0.0128,0.0129,0.0137,0.0112,0.0151,0.0224,0.0294,0.0309,0.0346,0.0364,0.0439,0.0415,0.0521,0.0580,0.0591,0.0539,0.0546,0.0472,0.0314,0.0263,0.0232,0.0193,0.0062,0.0060,-0.0037,-0.0110,-0.0168,-0.0078,-0.0119,-0.0097,-0.0012,0.0026,0.0062,0.0070,0.0166,0.0180,0.0227,0.0246,0.0257,0.0232,0.0291,0.0301,0.0286,0.0310,0.0322,0.0339,0.0352,0.0341,0.0335,0.0367,0.0375,0.0360,0.0355,0.0357,0.0385)
> 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.02179667
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 0.001223746
> (armose <- arm / armse)
[1] 17.81144
> (geo <- geomean(x))
[1] NaN
Warning message:
In log(x) : NaNs produced
> (har <- harmean(x))
[1] 0.02556666
> (qua <- quamean(x))
[1] 0.0255598
> (win <- winmean(x))
            [,1]         [,2]
 [1,] 0.02182833 0.0012112291
 [2,] 0.02178667 0.0011940861
 [3,] 0.02180167 0.0011826826
 [4,] 0.02180500 0.0011553171
 [5,] 0.02177167 0.0010770734
 [6,] 0.02173167 0.0010213951
 [7,] 0.02181333 0.0009568183
 [8,] 0.02184000 0.0008874943
 [9,] 0.02178000 0.0008736579
[10,] 0.02171333 0.0008637566
[11,] 0.02175917 0.0008488660
[12,] 0.02189917 0.0008177859
[13,] 0.02204000 0.0007906690
[14,] 0.02208667 0.0007787685
[15,] 0.02207417 0.0007703918
[16,] 0.02200750 0.0007574931
[17,] 0.02205000 0.0007344752
[18,] 0.02202000 0.0007303568
[19,] 0.02197250 0.0007199605
[20,] 0.02182250 0.0006839726
[21,] 0.02171750 0.0006707772
[22,] 0.02168083 0.0006662600
[23,] 0.02171917 0.0006618943
[24,] 0.02165917 0.0006545708
[25,] 0.02165917 0.0006496848
[26,] 0.02163750 0.0006470728
[27,] 0.02168250 0.0006315399
[28,] 0.02163583 0.0006045472
[29,] 0.02146667 0.0005849505
[30,] 0.02166667 0.0005465032
[31,] 0.02166667 0.0005408331
[32,] 0.02164000 0.0005378013
[33,] 0.02164000 0.0005258079
[34,] 0.02181000 0.0005022546
[35,] 0.02183917 0.0004993305
[36,] 0.02156917 0.0004627978
[37,] 0.02135333 0.0004267244
[38,] 0.02125833 0.0004102991
[39,] 0.02129083 0.0004070079
[40,] 0.02135750 0.0004003127
> (tri <- trimean(x))
            [,1]         [,2]
 [1,] 0.02180763 0.0011577014
 [2,] 0.02178621 0.0010971429
 [3,] 0.02178596 0.0010389846
 [4,] 0.02178036 0.0009776828
 [5,] 0.02177364 0.0009173071
 [6,] 0.02177407 0.0008723279
 [7,] 0.02178208 0.0008358669
 [8,] 0.02177692 0.0008095009
 [9,] 0.02176765 0.0007937647
[10,] 0.02176600 0.0007785667
[11,] 0.02177245 0.0007631115
[12,] 0.02177396 0.0007480622
[13,] 0.02176064 0.0007357900
[14,] 0.02173261 0.0007257270
[15,] 0.02169889 0.0007158392
[16,] 0.02166477 0.0007056201
[17,] 0.02163488 0.0006956082
[18,] 0.02160000 0.0006870138
[19,] 0.02156585 0.0006775383
[20,] 0.02153375 0.0006678464
[21,] 0.02151154 0.0006611350
[22,] 0.02149605 0.0006548080
[23,] 0.02148243 0.0006477780
[24,] 0.02146528 0.0006398667
[25,] 0.02145143 0.0006313210
[26,] 0.02143676 0.0006216672
[27,] 0.02142273 0.0006103833
[28,] 0.02140469 0.0005989095
[29,] 0.02138871 0.0005888220
[30,] 0.02138333 0.0005792951
[31,] 0.02136379 0.0005728940
[32,] 0.02134286 0.0005654952
[33,] 0.02132222 0.0005564856
[34,] 0.02130000 0.0005468682
[35,] 0.02126400 0.0005381431
[36,] 0.02122292 0.0005271212
[37,] 0.02119783 0.0005196012
[38,] 0.02118636 0.0005159260
[39,] 0.02118095 0.0005132307
[40,] 0.02117250 0.0005092225
> (midr <- midrange(x))
[1] 0.02115
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 0.02113548 0.02138333 0.02113548 0.02138333 0.02138333 0.02113548 0.02138333
[8] 0.02126667
> postscript(file="/var/www/rcomp/tmp/1asf41323889945.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/20uem1323889945.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/3enf91323889945.tab") 
> 
> try(system("convert tmp/1asf41323889945.ps tmp/1asf41323889945.png",intern=TRUE))
character(0)
> try(system("convert tmp/20uem1323889945.ps tmp/20uem1323889945.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.008   0.184   1.174