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(285 + ,NA + ,263.074978758798 + ,284.715000203273 + ,254.775 + ,215.375 + ,285 + ,251.150995958404 + ,273.669961695715 + ,224.125 + ,313.725 + ,278.0375 + ,240.148887802906 + ,251.121344746975 + ,288.75 + ,256.7625 + ,281.60625 + ,230.023039795889 + ,274.792708736568 + ,257.75 + ,273.79 + ,279.121875 + ,220.685960664563 + ,268.316598141974 + ,274.15 + ,173.0125 + ,278.5886875 + ,212.085468785266 + ,269.853519048767 + ,229.75 + ,174.875 + ,268.03106875 + ,204.143913113455 + ,235.891289968433 + ,214.05 + ,258.2625 + ,258.715461875 + ,196.82024483039 + ,214.685176171976 + ,92.75 + ,222.65 + ,258.6701656875 + ,190.082767981171 + ,229.822804324978 + ,206 + ,231.7375 + ,255.06814911875 + ,183.871806487569 + ,227.333358752735 + ,260.15 + ,150.778 + ,252.735084206875 + ,178.151865687297 + ,228.857278044667 + ,239.15 + ,144.1375 + ,242.539375786188 + ,172.868740338206 + ,201.763796135814 + ,178.4 + ,136.15 + ,232.699188207569 + ,167.995011999447 + ,181.769285845826 + ,222.875 + ,152.875 + ,223.044269386812 + ,163.498250062763 + ,165.941563187890 + ,160.75 + ,238.375 + ,216.027342448131 + ,159.353396166068 + ,161.40829991438 + ,160 + ,147.8 + ,218.262108203318 + ,155.548258305226 + ,188.11152745539 + ,130.025 + ,35.425 + ,211.215897382986 + ,152.031367572698 + ,174.125600162880 + ,137.25 + ,80.375 + ,193.636807644687 + ,148.767298136698 + ,126.004166756322 + ,77.3 + ,143.375 + ,182.310626880219 + ,145.755625456846 + ,110.173380876921 + ,105.05 + ,194.8875 + ,178.417064192197 + ,142.984891871717 + ,121.692499793419 + ,102.75 + ,190.43 + ,180.064107772977 + ,140.440637576800 + ,147.087096517147 + ,190.375 + ,122.525 + ,181.100696995679 + ,138.099808990679 + ,162.124674010416 + ,136.275 + ,153.125 + ,175.243127296111 + ,135.934098031119 + ,148.385789010579 + ,155.525 + ,79.6 + ,173.031314566500 + ,133.942182198235 + ,150.030031637519 + ,52.75 + ,182.8625 + ,163.688183109850 + ,132.093720050010 + ,125.594727004357 + ,131.625) + ,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.5074257 0 1.0000000 1.1140995 0.2774790 2 0.0100000 1 0.4925743 0.6852633 0.1686085 > postscript(file="/var/www/html/freestat/rcomp/tmp/1ip611274872811.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/2ip611274872811.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 285.0000 240.5190 44.4810 2 215.3750 240.5190 -25.1440 3 313.7250 240.5190 73.2060 4 256.7625 240.5190 16.2435 5 273.7900 240.5190 33.2710 6 173.0125 240.5190 -67.5065 7 174.8750 240.5190 -65.6440 8 258.2625 240.5190 17.7435 9 222.6500 240.5190 -17.8690 10 231.7375 240.5190 -8.7815 11 150.7780 143.5147 7.2633 12 144.1375 143.5147 0.6228 13 136.1500 143.5147 -7.3647 14 152.8750 143.5147 9.3603 15 238.3750 143.5147 94.8603 16 147.8000 143.5147 4.2853 17 35.4250 143.5147 -108.0897 18 80.3750 143.5147 -63.1397 19 143.3750 143.5147 -0.1397 20 194.8875 143.5147 51.3728 21 190.4300 143.5147 46.9153 22 122.5250 143.5147 -20.9897 23 153.1250 143.5147 9.6103 24 79.6000 143.5147 -63.9147 25 182.8625 143.5147 39.3478 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/freestat/rcomp/tmp/3sy5l1274872811.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/4wzma1274872811.tab") > > try(system("convert tmp/1ip611274872811.ps tmp/1ip611274872811.png",intern=TRUE)) character(0) > try(system("convert tmp/2ip611274872811.ps tmp/2ip611274872811.png",intern=TRUE)) character(0) > try(system("convert tmp/3sy5l1274872811.ps tmp/3sy5l1274872811.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.363 0.682 1.527