source('/home/pw/wessanet/cretab')



myrfcuid = ''

x <- c(20,25,15,15,25,25,25,21,30,25,20,40,13,30,25,20,25,20,25,20,20,15,15,12,20,5,20,15,25,22,20,22,25,20,20,35,30,25,20,20,20,25,25,15,20,35,25,25,30,23,10,22,25,25,22,30,20,25,25,22,25,25,25,22,25,12,18,20,20,22,30,25,22,20,50,30,25,20,30,22,25,30,22,25,22,22,25,25,25,20,22,15,20,30,20,25,30,35,22,12,30,15,10,30,9,25,20,20,35,25,35,30,12,25,15,25,25,20,20,6,15,40,20,40,25,25,20,15,15,22,24,22,20,25,25,25,35,40,20,22,22,20,25,25,18,25,20,25,30,20,22,35,22,25,25,25,25,22,23,35,15,25,18,22,25,25,28,30,20,25,25,30,22,30,10,10,25,20,22,25,25,15,22,25,25,28,22,30,25,20,25,25,20,30,20,30,50,19,20,28,20,25,35,25,25,15,16,20,20,25,30,20,25,25,25,20,20,25,25,30,22,20,25,25,18,18,20,25,25,30,25,20,25,20,20,20,22,18,22,20,15,25,25,20,25,15,22,25,25,15,12,25,30,22,15,22,25,12,18,30,25,25,40,24,25,15,25,20,25,25,25,20,30,20,25,30,22,25,25,25,50,19,50,25,35,20,20,20,20,20,25,25,25,20,20,20,20,25,18,25,22,22,30,30,8,20,25,30,50,22,20,10,25,25,25,25,18,25,20,25,30,18,20,25,22,22,20,20,25,20,20,20,20,25,20,10,20,25,30,25,50,30,30,50,15,25,25,22,20,22,30,25,18,22,22,30,40,25,20,10,20,9,15,20,15,20,30,12,15,12,20,15,12,25,20,25,25,25,30,20,25,15,15,22,10,15,10,20,25,20,20,38,20,20,20,40,25,25,30,25,10,20,25,12,15,25,20,22,22,20,25,25,25,15,40,20,20,16,25,15,20,25,20,30,50,20,25,20,30,30,25,25,12,25,25,25,20,20,20,15,20,25,15,25,50,30,20,20,25,12,15,20,20,35,22,15,18,30,22,12,12,20,20,15,25,15,20,20,25,18,30,20,25,25,25,20,20,25,20,22,15,15,22,20,10,25,20,20,15,12,20,5,20,15,15,25,25,25,15,25,22,25,20,18,22,25,35,25,25,25,35,30,22,30,50,15,25,24,20,25,25,25,12,15,22,25,25,25,25,15,20,20,15,35,30,20,22,65,20,25,22,20,25,25,20,25,15,20,12,15,10,25,15,30,35,25,25,25,25,25,40,40,25,25,20,25,25,22,25,30,25,25,30,25,25,30,25,25,20,22,22,20,25,22,25,22,40,25,25,25,22,20,35,20,35,25,22,25,25,25,25,25,40,25,30,25,20,25,25,30,22,22,20,15,15,25,25,20,20,15,25,15,20,22,25,15,15,18,5,15,25,18,40,25,25,20,30,20,25,25,25,22,22,25,25,30,25,25,25,25,20,20,25,25,25,25,20,30,25,22,30,20,20,30,25,25,30,20,25,25,24,25,30,18,15,22,22,25,22,22,25,15,20,22,18,35,20,20,20,25,25,30,15,25,22,26,25,20,25,25,25,22,25,25,20,22,30,15,30,25,20,25,25,35,22,20,25,20,20,18,20,22,25,10,20,25,20,20,30,25,20,15,20,25,10,20,25,22,22,25,25,15,25,20,10,25,16,25,35,25,15,25,25,30,25,10,22,20,25,20,20,25,22,18,30,19,25,20,25,20,25,20,22,12,30,12,22,25,25,25,25,30,30,10,22,22,25,20,22,20,25,20,15,25,20,25,20,30,15,40,25,20,22,22,30,20,40,20,25,20,25,20,50,50,25,25,40,30,22,30,20,25,25,30,25,25,20,18,18,28,25,22,15,40,40,12,12,18,12,25,26,18,25,22,15,25,15,15,15,25,15,12,22,20,20,25,20,12,9,15,12,15,25,20,20,15,15,30,21,25,22,22,50,15,25,15,25,22,18,50,20,50,20,20,30,25,20,22,25,50,40,25,25,25,25,30,40,25,30,20)
ylimmax = ''
ylimmin = ''
main = 'Robustness of Central Tendency'
#'GNU S' R Code compiled by R2WASP v. 1.2.327 (Thu, 20 Jul 2017 08:33:00 +0200)
#Author: root
#To cite this work: Wessa, P., (2017), Central Tendency (v1.0.7) in Free Statistics Software (v$_version), Office for Research Development and Education, URL https://www.wessa.net/rwasp_centraltendency.wasp/
#Source of accompanying publication: Office for Research, Development, and Education
#
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))
sqrtn <- sqrt(length(x))
(armse <- sd(x) / sqrtn)
(armose <- arm / armse)
(geo <- geomean(x))
(har <- harmean(x))
(qua <- quamean(x))
(win <- winmean(x))
(tri <- trimean(x))
(midr <- midrange(x))
midm <- array(NA,dim=8)
for (j in 1:8) midm[j] <- midmean(x,j)
midm
postscript(file="/home/pw/wessanet/rcomp/tmp/1ydj01508229913.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()
postscript(file="/home/pw/wessanet/rcomp/tmp/2yxcz1508229913.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()

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,'Arithmetic Mean',header=TRUE)
a<-table.element(a,signif(arm,6))
a<-table.element(a, signif(armse,6))
a<-table.element(a,signif(armose,6))
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Geometric Mean',header=TRUE)
a<-table.element(a,signif(geo,6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Harmonic Mean',header=TRUE)
a<-table.element(a,signif(har,6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Quadratic Mean',header=TRUE)
a<-table.element(a,signif(qua,6))
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, mylabel,header=TRUE)
a<-table.element(a,signif(win[j,1],6))
a<-table.element(a,signif(win[j,2],6))
a<-table.element(a,signif(win[j,1]/win[j,2],6))
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, mylabel,header=TRUE)
a<-table.element(a,signif(tri[j,1],6))
a<-table.element(a,signif(tri[j,2],6))
a<-table.element(a,signif(tri[j,1]/tri[j,2],6))
a<-table.row.end(a)
}
a<-table.row.start(a)
a<-table.element(a, 'Median',header=TRUE)
a<-table.element(a,signif(median(x),6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a, 'Midrange',header=TRUE)
a<-table.element(a,signif(midr,6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'Weighted Average at Xnp',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[1],6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'Weighted Average at X(n+1)p',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[2],6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'Empirical Distribution Function',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[3],6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'Empirical Distribution Function - Averaging',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[4],6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'Empirical Distribution Function - Interpolation',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[5],6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'Closest Observation',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[6],6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'True Basic - Statistics Graphics Toolkit',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[7],6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.row.start(a)
mymid <-  'Midmean'
mylabel <- paste(mymid,'MS Excel (old versions)',sep=' - ')
a<-table.element(a,mylabel,header=TRUE)
a<-table.element(a,signif(midm[8],6))
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,signif(length(x),6))
a<-table.element(a,'')
a<-table.element(a,'')
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file="/home/pw/wessanet/rcomp/tmp/3m1ak1508229913.tab") 

try(system("convert /home/pw/wessanet/rcomp/tmp/1ydj01508229913.ps /home/pw/wessanet/rcomp/tmp/1ydj01508229913.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/2yxcz1508229913.ps /home/pw/wessanet/rcomp/tmp/2yxcz1508229913.png",intern=TRUE))