R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(178600
+ ,181600
+ ,178500
+ ,176700
+ ,183500
+ ,175900
+ ,179800
+ ,186500
+ ,182700
+ ,182800
+ ,178600
+ ,183300
+ ,182900
+ ,191400
+ ,189300
+ ,192200
+ ,187900
+ ,193900
+ ,189100
+ ,193100
+ ,194800
+ ,200200
+ ,211500
+ ,202100
+ ,200300
+ ,199200
+ ,204900
+ ,207300
+ ,200000
+ ,197700
+ ,202200
+ ,200200
+ ,208300
+ ,215100
+ ,210700
+ ,208100
+ ,209000
+ ,211000
+ ,210200
+ ,205500
+ ,211400
+ ,211700
+ ,209300
+ ,207500
+ ,203300
+ ,207100
+ ,206900
+ ,228700
+ ,226900
+ ,265000
+ ,227100
+ ,228100
+ ,226500
+ ,225200
+ ,217800
+ ,221300
+ ,215300
+ ,231300
+ ,227100
+ ,237800
+ ,230200
+ ,233400
+ ,231100
+ ,237200
+ ,243700
+ ,239700
+ ,248400
+ ,241000
+ ,254500
+ ,242800
+ ,268300
+ ,253900
+ ,262100
+ ,264100
+ ,261000
+ ,269300
+ ,260400
+ ,263200
+ ,279200
+ ,272200
+ ,269200
+ ,289600
+ ,283200
+ ,284300
+ ,283000
+ ,289100
+ ,289600
+ ,289100
+ ,287400
+ ,279600
+ ,289300
+ ,295000
+ ,299600
+ ,293600
+ ,294400
+ ,290200
+ ,301000
+ ,307900
+ ,298800
+ ,310300
+ ,293900
+ ,305000
+ ,311300
+ ,317300
+ ,296200
+ ,306800
+ ,291800
+ ,301900
+ ,314600
+ ,321500
+ ,329400
+ ,311700
+ ,309700
+ ,306500
+ ,307100
+ ,301300
+ ,292200
+ ,310100
+ ,316800
+ ,284400
+ ,284600
+ ,301200
+ ,287600
+ ,314300
+ ,298200
+ ,299400
+ ,301900
+ ,265500
+ ,287100
+ ,274000
+ ,290100
+ ,263100
+ ,245200
+ ,258600
+ ,259800
+ ,269800
+ ,274600
+ ,274800
+ ,271100
+ ,257800
+ ,290300
+ ,262200
+ ,270000
+ ,290600)
> ylimmax = ''
> ylimmin = ''
> main = 'Robuustheid Huisprijzen USA'
> #'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] 249466.0
> sqrtn <- sqrt(length(x))
> (armse <- sd(x) / sqrtn)
[1] 3691.803
> (armose <- arm / armse)
[1] 67.57294
> (geo <- geomean(x))
[1] 245446.2
> (har <- harmean(x))
[1] 241369.3
> (qua <- quamean(x))
[1] 253342.2
> (win <- winmean(x))
          [,1]     [,2]
 [1,] 249416.7 3683.124
 [2,] 249383.3 3671.946
 [3,] 249375.0 3670.322
 [4,] 249313.9 3662.588
 [5,] 249345.1 3655.707
 [6,] 249311.8 3632.594
 [7,] 249345.8 3623.339
 [8,] 249295.8 3616.026
 [9,] 249289.6 3613.750
[10,] 249289.6 3606.928
[11,] 249167.4 3589.053
[12,] 249350.7 3549.974
[13,] 249450.0 3531.369
[14,] 249537.5 3513.933
[15,] 249402.1 3493.877
[16,] 249291.0 3428.412
[17,] 249385.4 3417.310
[18,] 249422.9 3396.168
[19,] 249515.3 3382.569
[20,] 249612.5 3365.297
[21,] 249831.2 3296.137
[22,] 250029.9 3267.772
[23,] 250061.8 3243.830
[24,] 249995.1 3229.762
[25,] 249647.9 3193.945
[26,] 249449.3 3170.060
[27,] 249674.3 3122.496
[28,] 249596.5 3110.713
[29,] 249757.6 3081.218
[30,] 249799.3 3017.457
[31,] 249842.4 2995.603
[32,] 249886.8 2937.598
[33,] 249863.9 2926.257
[34,] 249887.5 2919.155
[35,] 249911.8 2911.853
[36,] 249936.8 2884.565
[37,] 249988.2 2879.357
[38,] 250093.8 2853.091
[39,] 250120.8 2839.713
[40,] 250370.8 2814.672
[41,] 250086.1 2759.614
[42,] 250115.3 2745.346
[43,] 250145.1 2724.955
[44,] 249411.8 2650.211
[45,] 249411.8 2638.162
[46,] 250466.0 2527.743
[47,] 250172.2 2487.881
[48,] 250938.9 2400.598
> (tri <- trimean(x))
          [,1]     [,2]
 [1,] 249421.1 3664.421
 [2,] 249425.7 3643.768
 [3,] 249447.8 3627.168
 [4,] 249473.5 3609.311
 [5,] 249516.4 3591.731
 [6,] 249553.8 3573.813
 [7,] 249598.5 3558.655
 [8,] 249639.1 3543.294
 [9,] 249688.1 3527.178
[10,] 249739.5 3509.334
[11,] 249792.6 3490.218
[12,] 249860.8 3471.103
[13,] 249912.7 3454.576
[14,] 249956.9 3438.038
[15,] 249994.7 3421.206
[16,] 250045.5 3404.205
[17,] 250107.3 3391.853
[18,] 250163.9 3378.522
[19,] 250219.8 3365.079
[20,] 250271.2 3350.652
[21,] 250317.6 3335.442
[22,] 250351.0 3324.443
[23,] 250372.4 3313.861
[24,] 250392.7 3303.194
[25,] 250418.1 3291.389
[26,] 250466.3 3280.307
[27,] 250528.9 3268.774
[28,] 250580.7 3259.013
[29,] 250639.5 3247.604
[30,] 250691.7 3236.057
[31,] 250743.9 3227.572
[32,] 250796.2 3218.281
[33,] 250848.7 3211.568
[34,] 250905.3 3202.968
[35,] 250963.5 3191.720
[36,] 251023.6 3177.384
[37,] 251085.7 3161.581
[38,] 251148.5 3141.543
[39,] 251209.1 3118.940
[40,] 251271.9 3091.682
[41,] 251324.2 3060.376
[42,] 251396.7 3027.896
[43,] 251472.4 2988.524
[44,] 251551.8 2941.603
[45,] 251681.5 2893.043
[46,] 251821.2 2832.959
[47,] 251906.0 2776.017
[48,] 252016.7 2709.360
> (midr <- midrange(x))
[1] 252650
> midm <- array(NA,dim=8)
> for (j in 1:8) midm[j] <- midmean(x,j)
> midm
[1] 250427.4 251023.6 250427.4 251023.6 251023.6 250427.4 251023.6 250963.5
> postscript(file="/var/www/html/rcomp/tmp/1nx6x1268409392.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/2otk21268409392.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> 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/38szi1268409392.tab") 
> 
> try(system("convert tmp/1nx6x1268409392.ps tmp/1nx6x1268409392.png",intern=TRUE))
character(0)
> try(system("convert tmp/2otk21268409392.ps tmp/2otk21268409392.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.794   0.326   0.961