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 <- array(list(192 + ,NA + ,201.992383256927 + ,206.967782536146 + ,212 + ,212.25 + ,192 + ,197.374595777152 + ,201.340376307467 + ,219.9 + ,191.8 + ,194.025 + ,204.248977354664 + ,210.239313211865 + ,225.9 + ,163.7625 + ,193.8025 + ,198.495922223946 + ,201.252485572609 + ,204.05 + ,272.025 + ,190.7985 + ,182.444540071016 + ,188.931303179668 + ,256.65 + ,284.575 + ,198.92115 + ,223.842424405143 + ,236.507656592625 + ,260.375 + ,301.6635 + ,207.486535 + ,251.908814573073 + ,242.02280020498 + ,221.75 + ,287.5375 + ,216.9042315 + ,274.901984222098 + ,249.532404318106 + ,123.05 + ,220.4375 + ,223.96755835 + ,280.741244492589 + ,243.324681715068 + ,277.375 + ,178.3 + ,223.614552515 + ,252.873030359752 + ,213.837340170125 + ,244.4 + ,284.8875 + ,219.0830972635 + ,218.410540574304 + ,195.319860969677 + ,223.8 + ,283.9875 + ,225.66353753715 + ,249.131587099083 + ,242.160129478196 + ,310.9 + ,238 + ,231.495933783435 + ,265.239575959338 + ,241.764621171334 + ,254.4 + ,216.275 + ,232.146340405092 + ,252.651330534283 + ,221.555245324865 + ,254.625 + ,162.875 + ,230.559206364582 + ,235.840709946355 + ,212.008114250888 + ,122.65 + ,185.95 + ,223.790785728124 + ,202.121012334719 + ,188.541288043735 + ,117.775 + ,193.7875 + ,220.006707155312 + ,194.647890420306 + ,198.681681578005 + ,187.15 + ,128.3275 + ,217.384786439781 + ,194.250277557744 + ,202.125899750263 + ,195.65 + ,83.925 + ,208.479057795803 + ,163.7853355103 + ,173.359262231158 + ,115.7 + ,177.15 + ,196.023652016222 + ,126.879419477262 + ,153.846420458440 + ,144.375 + ,142.3 + ,194.1362868146 + ,150.111000109239 + ,194.814489244242 + ,125.75 + ,120.5375 + ,188.95265813314 + ,146.501296840105 + ,179.499528695192 + ,115.65 + ,269.25 + ,182.111142319826 + ,134.502628164888 + ,169.935918108429 + ,193 + ,167.625 + ,190.825028087843 + ,196.773530970922 + ,235.288172646467 + ,184.35 + ,243.275 + ,188.505025279059 + ,183.303098751657 + ,190.628692996618 + ,194.325) + ,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.2987862 0 1.0000000 1.007148 0.2169807 2 0.0100000 1 0.7012138 1.244799 0.2464910 > postscript(file="/var/www/html/rcomp/tmp/16ann1274873413.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/rcomp/tmp/2hj481274873413.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 192.0000 182.0368 9.963167 2 212.2500 182.0368 30.213167 3 191.8000 247.1951 -55.395100 4 163.7625 182.0368 -18.274333 5 272.0250 247.1951 24.829900 6 284.5750 247.1951 37.379900 7 301.6635 247.1951 54.468400 8 287.5375 182.0368 105.500667 9 220.4375 247.1951 -26.757600 10 178.3000 247.1951 -68.895100 11 284.8875 247.1951 37.692400 12 283.9875 247.1951 36.792400 13 238.0000 247.1951 -9.195100 14 216.2750 247.1951 -30.920100 15 162.8750 182.0368 -19.161833 16 185.9500 182.0368 3.913167 17 193.7875 182.0368 11.750667 18 128.3275 182.0368 -53.709333 19 83.9250 182.0368 -98.111833 20 177.1500 182.0368 -4.886833 21 142.3000 182.0368 -39.736833 22 120.5375 182.0368 -61.499333 23 269.2500 182.0368 87.213167 24 167.6250 182.0368 -14.411833 25 243.2750 182.0368 61.238167 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/rcomp/tmp/3hj481274873413.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/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,'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/rcomp/tmp/4vb2z1274873413.tab") > > try(system("convert tmp/16ann1274873413.ps tmp/16ann1274873413.png",intern=TRUE)) character(0) > try(system("convert tmp/2hj481274873413.ps tmp/2hj481274873413.png",intern=TRUE)) character(0) > try(system("convert tmp/3hj481274873413.ps tmp/3hj481274873413.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.931 0.474 3.459