R version 2.6.1 (2007-11-26)
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(0.112099938198320
+ ,-7.5234317022524
+ ,-0.174508990247777
+ ,-2.7838976248972
+ ,3.32086811233621
+ ,-2.45317788349593
+ ,3.18406012068406
+ ,-5.44042146983815
+ ,0.649481017642372
+ ,-1.07642801808016
+ ,-6.95667668647274
+ ,-1.51981598559022
+ ,4.16227230153779
+ ,1.31150417256041
+ ,-3.04705830224746
+ ,-1.92422722261651
+ ,2.75950343482618
+ ,-6.41361862382443
+ ,-1.11530752560768
+ ,-1.56757720439110
+ ,-1.87571926492828
+ ,-0.248344227067880
+ ,5.19714930126309
+ ,0.523454484546629
+ ,1.7483538545492
+ ,-0.151122707042641
+ ,6.41455842806672
+ ,-2.33399246036476
+ ,-0.345905183770059
+ ,-1.93508355264285
+ ,5.12010225384436
+ ,-1.81633381269994
+ ,11.8083712468710
+ ,0.951312732337543
+ ,-2.00735217132620
+ ,1.39325355869821
+ ,5.89468248583685
+ ,2.47323836463487
+ ,-10.8055792141131
+ ,1.87601361620464
+ ,-4.67348222231403
+ ,-7.31601418786265
+ ,4.05419945496961
+ ,-0.73915923860179
+ ,4.19522631703279
+ ,-6.8416643186816
+ ,0.922135364442534
+ ,3.04248996617448
+ ,6.7374333015147
+ ,-1.80882009330620
+ ,1.70018100638115
+ ,-3.87944179177855
+ ,2.02934054306604
+ ,2.46312058350850
+ ,0.0958696906479143
+ ,2.06857942650907
+ ,0.451957173332116
+ ,-0.169740018337748
+ ,8.36055226069874
+ ,6.95713440198822
+ ,-0.737722119028973
+ ,10.8460464500845
+ ,10.4303191959713
+ ,-9.3254800323235
+ ,2.98350849441621
+ ,-3.37677921920559
+ ,-4.70835338252925
+ ,-1.60721490048834
+ ,-3.14479871850406
+ ,5.22904587501341
+ ,1.71870765649422
+ ,7.8616002868662
+ ,16.6634543338289
+ ,-1.86414944852601
+ ,-2.56744340647083
+ ,0.579191901975491
+ ,-3.42163983030900
+ ,-6.11260379447978
+ ,6.30214100944247
+ ,-3.35893281479719
+ ,-2.22878328457179
+ ,-1.44249975349172
+ ,8.89655891765128
+ ,4.94508794854377
+ ,-2.25423698595273
+ ,1.09067021401188
+ ,-3.32773833069336
+ ,0.423899895235991
+ ,-0.782494091513541
+ ,5.71247870563437
+ ,3.12159570925235
+ ,-4.26681435448836
+ ,6.29041863604803)
> 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] 0.5548563
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 0.4948582
> (armose <- arm / armse)
[1] 1.121243
> (geo <- geomean(x))
[1] NaN
Warning message:
In log(x) : NaNs produced
> (har <- harmean(x))
[1] -114.2745
> (qua <- quamean(x))
[1] 4.778833
> (win <- winmean(x))
            [,1]      [,2]
 [1,] 0.51856617 0.4750466
 [2,] 0.53662474 0.4618103
 [3,] 0.52990507 0.4573467
 [4,] 0.47939227 0.4398457
 [5,] 0.45675817 0.4328529
 [6,] 0.45218357 0.4216063
 [7,] 0.40676252 0.4051241
 [8,] 0.44568564 0.3920339
 [9,] 0.48528498 0.3755102
[10,] 0.47694667 0.3728832
[11,] 0.52366065 0.3656088
[12,] 0.52258148 0.3499902
[13,] 0.56110596 0.3372647
[14,] 0.49508433 0.3245786
[15,] 0.49281818 0.3233912
[16,] 0.48492957 0.3206074
[17,] 0.48637828 0.3113215
[18,] 0.36016127 0.2870366
[19,] 0.40719274 0.2792464
[20,] 0.43050056 0.2701537
[21,] 0.29071150 0.2437138
[22,] 0.28654272 0.2359390
[23,] 0.29081901 0.2315121
[24,] 0.27697332 0.2280251
[25,] 0.32064258 0.2189349
[26,] 0.27822164 0.2084070
[27,] 0.19826426 0.1974547
[28,] 0.20982259 0.1953725
[29,] 0.09040142 0.1798703
[30,] 0.09316812 0.1765416
[31,] 0.04456371 0.1702101
> (tri <- trimean(x))
            [,1]      [,2]
 [1,] 0.50267871 0.4564607
 [2,] 0.48607721 0.4349989
 [3,] 0.45906042 0.4185070
 [4,] 0.43322296 0.4014196
 [5,] 0.42029000 0.3879372
 [6,] 0.41191582 0.3744299
 [7,] 0.40401519 0.3617528
 [8,] 0.40354116 0.3510021
 [9,] 0.39700877 0.3413560
[10,] 0.38451305 0.3334534
[11,] 0.37240555 0.3246446
[12,] 0.35387231 0.3155665
[13,] 0.35387231 0.3076709
[14,] 0.30940169 0.3004486
[15,] 0.28982291 0.2940377
[16,] 0.26919060 0.2864348
[17,] 0.24793666 0.2776437
[18,] 0.22505217 0.2685386
[19,] 0.21236010 0.2618934
[20,] 0.19436662 0.2549298
[21,] 0.17283676 0.2477829
[22,] 0.16218336 0.2437563
[23,] 0.15099823 0.2398549
[24,] 0.13843463 0.2354492
[25,] 0.12595004 0.2302196
[26,] 0.12595004 0.2249533
[27,] 0.09269937 0.2200735
[28,] 0.08287201 0.2158130
[29,] 0.07082466 0.2101050
[30,] 0.06892221 0.2059627
[31,] 0.06649762 0.2005331
> (midr <- midrange(x))
[1] 2.928938
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 0.08642003 0.15099823 0.15099823 0.15099823 0.15099823 0.09922759 0.15099823
[8] 0.15099823
> postscript(file="/var/www/html/rcomp/tmp/10qs91197838198.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/2rffl1197838198.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/3g79w1197838199.tab") 
> 
> system("convert tmp/10qs91197838198.ps tmp/10qs91197838198.png")
> system("convert tmp/2rffl1197838198.ps tmp/2rffl1197838198.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.670   0.530   1.774