R version 2.8.0 (2008-10-20) Copyright (C) 2008 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(15.22,15.27,15.31,15.33,15.42,15.49,15.65,15.67,15.69,15.83,15.92,15.99,15.94,15.96,16.03,16.09,16.04,16.23,16.2,16.2,16.26,16.28,16.27,16.29,16.3,16.37,16.39,16.42,16.43,16.37,16.37,16.39,16.48,16.51,16.5,16.54,16.52,16.56,16.61,16.75,16.75,16.79,16.82,16.84,17.14,17.25,17.28,17.3,17.34,17.44,17.48,17.55,17.59,17.66,17.67,17.64,17.68,17.72,17.78,17.83,17.88,18.11,18.16,18.27,18.29,18.35,18.35,18.38,18.41,18.41,18.42,18.43,18.48,18.54,18.65,18.66,18.69,18.72,18.72,18.73,18.84,18.83,18.91,18.91,18.94,18.97,19,19.08,19.18,19.24,19.23,19.25,19.3,19.33,19.35,19.35,19.31,19.47,19.7,19.76,19.9,19.97,20.1,20.26,20.44,20.43,20.57,20.6,20.69,20.93,20.98,21.11,21.14,21.16,21.32,21.32,21.48,21.58,21.74,21.75,21.81,21.89,22.21,22.37,22.47,22.51,22.55,22.61,22.58,22.85,22.93,22.98) > ylimmax = '' > ylimmin = '' > main = 'Centrummaten gem. pr. gebak. tong of for' > #'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] 18.42023 > sqrtn <- sqrt(length(x)) > (armse <- sd(x) / sqrtn) [1] 0.1841135 > (armose <- arm / armse) [1] 100.0482 > (geo <- geomean(x)) [1] 18.30289 > (har <- harmean(x)) [1] 18.18902 > (qua <- quamean(x)) [1] 18.54037 > (win <- winmean(x)) [,1] [,2] [1,] 18.42023 0.18399238 [2,] 18.41962 0.18368864 [3,] 18.41462 0.18264957 [4,] 18.41644 0.18214351 [5,] 18.41795 0.18161603 [6,] 18.42341 0.18042669 [7,] 18.42235 0.17993705 [8,] 18.41750 0.17876386 [9,] 18.41614 0.17585292 [10,] 18.39871 0.17120886 [11,] 18.39371 0.16999467 [12,] 18.39008 0.16896230 [13,] 18.39205 0.16848990 [14,] 18.37932 0.16548917 [15,] 18.36909 0.16370273 [16,] 18.35576 0.16027042 [17,] 18.36992 0.15876688 [18,] 18.34811 0.15571532 [19,] 18.34955 0.15486598 [20,] 18.34955 0.15376871 [21,] 18.33045 0.15079889 [22,] 18.32379 0.14951172 [23,] 18.28371 0.14387476 [24,] 18.26917 0.14160938 [25,] 18.27674 0.13950055 [26,] 18.25114 0.13632425 [27,] 18.24909 0.13607377 [28,] 18.21727 0.13127602 [29,] 18.18212 0.12716381 [30,] 18.15939 0.12306874 [31,] 18.14530 0.12098167 [32,] 18.12348 0.11596202 [33,] 18.11348 0.11381489 [34,] 18.05682 0.10740172 [35,] 18.02765 0.10396577 [36,] 18.03311 0.10336313 [37,] 18.03311 0.10220149 [38,] 18.04174 0.10006611 [39,] 18.08015 0.09531544 [40,] 18.06500 0.09384502 [41,] 18.07432 0.09222054 [42,] 18.08068 0.09090080 [43,] 18.07091 0.08864672 [44,] 18.13758 0.07506726 > (tri <- trimean(x)) [,1] [,2] [1,] 18.40977 0.18193836 [2,] 18.39898 0.17967733 [3,] 18.38817 0.17735977 [4,] 18.37879 0.17521717 [5,] 18.36861 0.17300457 [6,] 18.35775 0.17069037 [7,] 18.34551 0.16838745 [8,] 18.33302 0.16592660 [9,] 18.32079 0.16340677 [10,] 18.30830 0.16109363 [11,] 18.29745 0.15925161 [12,] 18.28676 0.15736543 [13,] 18.27604 0.15539167 [14,] 18.26471 0.15323040 [15,] 18.25412 0.15119933 [16,] 18.24400 0.14913813 [17,] 18.23459 0.14724035 [18,] 18.22365 0.14525259 [19,] 18.21394 0.14336961 [20,] 18.20370 0.14131276 [21,] 18.19300 0.13908390 [22,] 18.18318 0.13689453 [23,] 18.17337 0.13453010 [24,] 18.16583 0.13250621 [25,] 18.15890 0.13043554 [26,] 18.15113 0.12827038 [27,] 18.14462 0.12614454 [28,] 18.13789 0.12368215 [29,] 18.13284 0.12141498 [30,] 18.12972 0.11927516 [31,] 18.12786 0.11726364 [32,] 18.12676 0.11513332 [33,] 18.12697 0.11324170 [34,] 18.12781 0.11123862 [35,] 18.13226 0.10965077 [36,] 18.13883 0.10812634 [37,] 18.14552 0.10630663 [38,] 18.15268 0.10420026 [39,] 18.15981 0.10190310 [40,] 18.16500 0.09983149 [41,] 18.17160 0.09741364 [42,] 18.17813 0.09461740 [43,] 18.18478 0.09127055 [44,] 18.19273 0.08736836 > (midr <- midrange(x)) [1] 19.1 > midm <- array(NA,dim=8) > for (j in 1:8) midm[j] <- midmean(x,j) > midm [1] 18.10239 18.12697 18.10239 18.12697 18.12697 18.10239 18.12697 18.12676 > postscript(file="/var/www/html/freestat/rcomp/tmp/1kyhk1250534077.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/freestat/rcomp/tmp/2dinw1250534077.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/370lr1250534077.tab") > > system("convert tmp/1kyhk1250534077.ps tmp/1kyhk1250534077.png") > system("convert tmp/2dinw1250534077.ps tmp/2dinw1250534077.png") > > > proc.time() user system elapsed 1.099 0.519 1.197