R version 2.6.0 (2007-10-03)
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(622.755354196211
+ ,419.557064448606
+ ,337.667064448595
+ ,587.237064448597
+ ,-78.5029355514001
+ ,568.837064448604
+ ,111.996766841012
+ ,336.956766841012
+ ,-410.723233158992
+ ,506.376766841009
+ ,-207.193233158988
+ ,342.366766841013
+ ,395.671231421387
+ ,256.272941673806
+ ,221.182941673811
+ ,-171.54705832619
+ ,-252.987058326191
+ ,-506.447058326191
+ ,-1043.78735593379
+ ,-941.627355933786
+ ,-841.907355933784
+ ,-600.707355933784
+ ,-876.677355933786
+ ,-265.917355933785
+ ,-632.912891353406
+ ,-85.1111811009849
+ ,-339.801181100985
+ ,-361.331181100985
+ ,-139.271181100985
+ ,302.068818899015
+ ,-360.67147870858
+ ,645.48852129142
+ ,410.608521291422
+ ,-300.891478708579
+ ,1072.53852129142
+ ,-225.101478708580
+ ,109.802985871800
+ ,706.304696124221
+ ,382.814696124222
+ ,358.684696124221
+ ,1231.64469612422
+ ,593.68469612422
+ ,1026.44439851663
+ ,641.904398516626
+ ,79.9243985166273
+ ,249.724398516626
+ ,720.354398516625
+ ,-9.28560148337375
+ ,-389.081136902994
+ ,-315.179426650573
+ ,-800.369426650573
+ ,-673.499426650573
+ ,-920.739426650575
+ ,-1487.59942665057
+ ,370.863251817773
+ ,198.223251817774
+ ,-9.05674818222588
+ ,870.343251817775
+ ,572.773251817774
+ ,-117.366748182226
+ ,593.537716398154
+ ,-31.160573349427
+ ,87.7494266505743
+ ,798.919426650575
+ ,-181.820573349426
+ ,199.819426650573
+ ,314.47912904298
+ ,58.83912904298
+ ,-547.640870957020
+ ,49.7591290429804
+ ,-771.71087095702
+ ,-701.150870957021
+ ,-141.946406376641
+ ,-1102.04469612422
+ ,-506.434696124219
+ ,209.935303875781
+ ,-826.90469612422
+ ,-268.264696124222
+ ,-677.204993731815
+ ,-493.544993731816
+ ,-136.324993731815
+ ,-191.124993731815
+ ,-1176.29499373181
+ ,438.965006268185
+ ,-431.030529151435
+ ,-146.428818899016
+ ,86.0811811009855
+ ,-81.7488188990137
+ ,-65.288818899014
+ ,493.651181100984
+ ,-478.489116506611
+ ,-296.029116506610
+ ,53.1908834933911
+ ,305.190883493390
+ ,-424.979116506609
+ ,84.5808834933917
+ ,-889.51465192623
+ ,-80.7129416738088
+ ,432.597058326191
+ ,-769.932941673809
+ ,828.527058326191
+ ,1028.86705832619
+ ,638.626760718596
+ ,307.086760718597
+ ,1529.50676071860
+ ,-202.193239281403
+ ,800.236760718596
+ ,299.296760718597
+ ,-377.998774701024
+ ,378.502935551396
+ ,98.5129355513976
+ ,103.282935551397
+ ,405.342935551397
+ ,-924.617064448604
+ ,97.7426379438026
+ ,-457.297362056198
+ ,-127.577362056198
+ ,-686.477362056198
+ ,290.952637943803
+ ,153.612637943802
+ ,1140.71710252418)
> 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] 6.796767e-14
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 51.40381
> (armose <- arm / armse)
[1] 1.322230e-15
> (geo <- geomean(x))
[1] NaN
Warning message:
In log(x) : NaNs produced
> (har <- harmean(x))
[1] -611.2891
> (qua <- quamean(x))
[1] 563.1005
> (win <- winmean(x))
             [,1]     [,2]
 [1,]   0.1110940 50.28515
 [2,]  -0.1645638 49.75523
 [3,]  -0.4105450 49.17902
 [4,]   1.5229603 48.34667
 [5,]   2.1257550 48.21552
 [6,]  -5.4225181 46.89584
 [7,]  -6.0352448 46.23487
 [8,]  -7.0569309 45.82218
 [9,]  -4.5687161 45.40451
[10,]  -9.8218044 44.28873
[11,]  -8.6867528 43.74665
[12,] -11.8759382 42.52112
[13,] -12.0699920 42.44316
[14,]  -4.4909667 41.24149
[15,]  -4.6394664 40.73245
[16,]  -7.2574221 40.07459
[17,]  -6.7574553 39.99959
[18,]  -1.6570759 39.05600
[19,]   1.1288144 38.09778
[20,]   9.2495248 36.88647
[21,]   5.5586472 34.65742
[22,]   3.2471520 34.38031
[23,]  -4.6976320 32.87409
[24,]  -2.9744064 32.36437
[25,]  -1.2901583 31.53150
[26,]   2.4311271 30.64500
[27,]   2.6064769 30.35601
[28,]   3.6672787 29.72403
[29,]   5.7729090 28.78874
[30,]   7.4515706 28.35948
[31,]   9.7645062 27.66580
[32,]   6.7181978 27.29807
[33,]   7.9597528 26.18794
[34,]  13.5576848 25.28949
[35,]  17.4851043 24.82276
[36,]  12.2441959 23.95761
[37,]  18.4736664 22.81974
[38,]  18.6154481 22.67989
[39,]  21.7767794 22.14044
[40,]  30.0787700 21.10196
> (tri <- trimean(x))
            [,1]     [,2]
 [1,] -0.3521625 49.07243
 [2,] -0.8312567 47.72600
 [3,] -1.1819952 46.54180
 [4,] -1.4573505 45.46273
 [5,] -2.2695523 44.53106
 [6,] -3.2453912 43.52730
 [7,] -2.8350605 42.71942
 [8,] -2.3082274 41.95451
 [9,] -1.6109056 41.17599
[10,] -1.2171817 40.37640
[11,] -0.1655056 39.66338
[12,]  0.8008214 38.94483
[13,]  2.1463371 38.31697
[14,]  3.5691459 37.61433
[15,]  4.3346668 36.98888
[16,]  5.1480527 36.34482
[17,]  6.2264021 35.69594
[18,]  7.3136317 34.95997
[19,]  8.0401750 34.24227
[20,]  8.5835633 33.54272
[21,]  8.5325624 32.89422
[22,]  8.7551003 32.42956
[23,]  9.1590165 31.91728
[24,] 10.1576195 31.50324
[25,] 11.0901167 31.07415
[26,] 11.9585302 30.66735
[27,] 12.6203073 30.29170
[28,] 13.3107196 29.87078
[29,] 13.9722027 29.44174
[30,] 14.5330363 29.04334
[31,] 14.5330363 28.60837
[32,] 15.3768242 28.16676
[33,] 15.9721048 27.66523
[34,] 16.5264185 27.20674
[35,] 16.7335792 26.76585
[36,] 16.6805562 26.27221
[37,] 16.9978136 25.78053
[38,] 16.8905595 25.34027
[39,] 16.7628290 24.77270
[40,] 16.3834119 24.12047
> (midr <- midrange(x))
[1] 20.95367
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1]  8.466871 14.533036 14.533036 14.533036 14.533036  8.023130 14.533036
[8] 14.533036
> postscript(file="/var/www/html/rcomp/tmp/1hexf1198253050.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/2gisg1198253050.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/32wy81198253051.tab") 
> 
> system("convert tmp/1hexf1198253050.ps tmp/1hexf1198253050.png")
> system("convert tmp/2gisg1198253050.ps tmp/2gisg1198253050.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.118   0.329   1.223