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 <- array(list(341.25 + ,NA + ,333.775928269088 + ,340.90875028286 + ,381.75 + ,303.6875 + ,341.25 + ,338.492504687478 + ,332.86089603337 + ,417.75 + ,357.5 + ,337.49375 + ,316.528510995532 + ,325.788548906474 + ,359.25 + ,295.075 + ,339.494375 + ,342.383916967292 + ,327.541311770725 + ,336 + ,386.5755 + ,335.0524375 + ,312.529223276755 + ,332.035063343418 + ,361.8 + ,455.6625 + ,340.20474375 + ,359.256754833345 + ,355.979142780996 + ,420.275 + ,424.926 + ,351.750519375 + ,420.094418776242 + ,422.853279739822 + ,400.9 + ,506.751 + ,359.0680674375 + ,423.143429126151 + ,447.538297256219 + ,248.15 + ,433.9 + ,373.83636069375 + ,475.904696528962 + ,463.452634199947 + ,480.125 + ,466.3375 + ,379.842724624375 + ,449.397276177252 + ,469.038817658043 + ,441.4 + ,496.7 + ,388.492202161938 + ,460.087548237949 + ,454.838864911353 + ,483 + ,464.45 + ,399.312981945744 + ,483.19214675098 + ,483.073307566238 + ,477.75 + ,385.375 + ,405.826683751169 + ,471.364755928644 + ,464.126997390433 + ,486.25 + ,381.875 + ,403.781515376053 + ,417.100190422978 + ,429.475306145968 + ,433.875 + ,219.6375 + ,401.590863838447 + ,394.871034955396 + ,366.780195970248 + ,286.4 + ,268.975 + ,383.395527454603 + ,284.288425042335 + ,314.305881579131 + ,231.875 + ,292.2875 + ,371.953474709142 + ,274.624758163045 + ,216.573154866353 + ,308.65 + ,181.025 + ,363.986877238228 + ,285.770981350566 + ,228.421842054758 + ,253.025 + ,277.625 + ,345.690689514405 + ,219.670140557742 + ,233.11741103145 + ,201.025 + ,166.75 + ,338.884120562965 + ,256.243046627923 + ,200.091184304044 + ,189.775 + ,266 + ,321.670708506668 + ,199.767699784603 + ,235.499484918374 + ,237.375 + ,189.25 + ,316.103637656002 + ,241.56415512763 + ,219.691637200487 + ,171.125 + ,226.35 + ,303.418273890401 + ,208.550863403060 + ,204.478173107192 + ,212.4 + ,158.75 + ,295.711446501361 + ,219.783159661820 + ,216.921637949675 + ,143.1 + ,218.8125 + ,282.015301851225 + ,181.267666104251 + ,235.210210182817 + ,205.125) + ,dim=c(5 + ,25) + ,dimnames=list(c('Actuals' + ,'Croston' + ,'ETS' + ,'Arima' + ,'Kaneka') + ,1:25)) > y <- array(NA,dim=c(5,25),dimnames=list(c('Actuals','Croston','ETS','Arima','Kaneka'),1:25)) > for (i in 1:dim(x)[1]) + { + for (j in 1:dim(x)[2]) + { + y[i,j] <- as.numeric(x[i,j]) + } + } > par2 = 'No' > par1 = '1' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Dr. Ian E. Holliday > #To cite this work: Ian E. Holliday, 2009, 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: > #Technical description: > library(rpart) > library(partykit) Loading required package: grid Loading required package: mvtnorm > par1 <- as.numeric(par1) > autoprune <- function ( tree, method='Minimum CV'){ + xerr <- tree$cptable[,'xerror'] + cpmin.id <- which.min(xerr) + if (method == 'Minimum CV Error plus 1 SD'){ + xstd <- tree$cptable[,'xstd'] + errt <- xerr[cpmin.id] + xstd[cpmin.id] + cpSE1.min <- which.min( errt < xerr ) + mycp <- (tree$cptable[,'CP'])[cpSE1.min] + } + if (method == 'Minimum CV') { + mycp <- (tree$cptable[,'CP'])[cpmin.id] + } + return (mycp) + } > conf.multi.mat <- function(true, new) + { + if ( all( is.na(match( levels(true),levels(new) ) )) ) + stop ( 'conflict of vector levels') + multi.t <- list() + for (mylev in levels(true) ) { + true.tmp <- true + new.tmp <- new + left.lev <- levels (true.tmp)[- match(mylev,levels(true) ) ] + levels(true.tmp) <- list ( mylev = mylev, all = left.lev ) + levels(new.tmp) <- list ( mylev = mylev, all = left.lev ) + curr.t <- conf.mat ( true.tmp , new.tmp ) + multi.t[[mylev]] <- curr.t + multi.t[[mylev]]$precision <- + round( curr.t$conf[1,1] / sum( curr.t$conf[1,] ), 2 ) + } + return (multi.t) + } > x <- t(y) > k <- length(x[1,]) > n <- length(x[,1]) > x1 <- cbind(x[,par1], x[,1:k!=par1]) > mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1]) > colnames(x1) <- mycolnames #colnames(x)[par1] > m <- rpart(as.data.frame(x1)) > par2 [1] "No" > if (par2 != 'No') { + mincp <- autoprune(m,method=par2) + print(mincp) + m <- prune(m,cp=mincp) + } > m$cptable CP nsplit rel error xerror xstd 1 0.6025141 0 1.0000000 1.130355 0.2000363 2 0.0100000 1 0.3974859 1.009521 0.2845941 > postscript(file="/var/www/html/freestat/rcomp/tmp/1sy5l1274873201.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(as.party(m),tp_args=list(id=FALSE)) > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/2sy5l1274873201.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plotcp(m) > dev.off() null device 1 > cbind(y=m$y,pred=predict(m),res=residuals(m)) y pred res 1 341.2500 394.6468 -53.396833 2 303.6875 394.6468 -90.959333 3 357.5000 394.6468 -37.146833 4 295.0750 394.6468 -99.571833 5 386.5755 394.6468 -8.071333 6 455.6625 394.6468 61.015667 7 424.9260 394.6468 30.279167 8 506.7510 394.6468 112.104167 9 433.9000 394.6468 39.253167 10 466.3375 394.6468 71.690667 11 496.7000 394.6468 102.053167 12 464.4500 394.6468 69.803167 13 385.3750 394.6468 -9.271833 14 381.8750 394.6468 -12.771833 15 219.6375 394.6468 -175.009333 16 268.9750 224.5825 44.392500 17 292.2875 224.5825 67.705000 18 181.0250 224.5825 -43.557500 19 277.6250 224.5825 53.042500 20 166.7500 224.5825 -57.832500 21 266.0000 224.5825 41.417500 22 189.2500 224.5825 -35.332500 23 226.3500 224.5825 1.767500 24 158.7500 224.5825 -65.832500 25 218.8125 224.5825 -5.770000 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/freestat/rcomp/tmp/3sy5l1274873201.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow=c(2,2)) > plot(myr,ylab='residuals') > plot(density(myr),main='Residual Kernel Density') > plot(myp,myr,xlab='predicted',ylab='residuals',main='Predicted vs Residuals') > plot(density(myp),main='Prediction Kernel Density') > par(op) > 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,'Model Performance',6,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'#',header=TRUE) > a<-table.element(a,'Complexity',header=TRUE) > a<-table.element(a,'split',header=TRUE) > a<-table.element(a,'relative error',header=TRUE) > a<-table.element(a,'CV error',header=TRUE) > a<-table.element(a,'CV S.D.',header=TRUE) > a<-table.row.end(a) > for (i in 1:length(m$cptable[,1])) { + a<-table.row.start(a) + a<-table.element(a,i,header=TRUE) + a<-table.element(a,round(m$cptable[i,'CP'],3)) + a<-table.element(a,m$cptable[i,'nsplit']) + a<-table.element(a,round(m$cptable[i,'rel error'],3)) + a<-table.element(a,round(m$cptable[i,'xerror'],3)) + a<-table.element(a,round(m$cptable[i,'xstd'],3)) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/freestat/rcomp/tmp/4p83c1274873201.tab") > > try(system("convert tmp/1sy5l1274873201.ps tmp/1sy5l1274873201.png",intern=TRUE)) character(0) > try(system("convert tmp/2sy5l1274873201.ps tmp/2sy5l1274873201.png",intern=TRUE)) character(0) > try(system("convert tmp/3sy5l1274873201.ps tmp/3sy5l1274873201.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.368 0.674 1.538