R version 2.10.1 (2009-12-14) 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(46 + ,NA + ,34.4268537070987 + ,45.95244183468 + ,31 + ,40.5 + ,46 + ,32.4378603466732 + ,40.45638364361 + ,50 + ,22.5 + ,45.45 + ,30.492240334061 + ,22.4728254462900 + ,20 + ,25 + ,43.155 + ,28.5853543104307 + ,24.9687672592201 + ,30 + ,22.25 + ,41.3395 + ,26.7195658658305 + ,22.2199590695251 + ,29 + ,7 + ,39.43055 + ,24.8931981461802 + ,6.98365087358009 + ,10.25 + ,11 + ,36.187495 + ,23.1020415201185 + ,10.9780926872601 + ,11 + ,50.25 + ,33.6687455 + ,21.3485568619093 + ,50.1872845185651 + ,15.25 + ,16.25 + ,35.32687095 + ,19.6418079270611 + ,16.2197263132451 + ,21 + ,32.5 + ,33.419183855 + ,17.9671794656517 + ,32.4519181330501 + ,15 + ,5.7525 + ,33.3272654695 + ,16.3327560407159 + ,5.72960743135642 + ,30 + ,7.75 + ,30.56978892255 + ,14.7288766581764 + ,7.72355174403519 + ,5.6 + ,14 + ,28.287810030295 + ,13.1590323347601 + ,27.3017308111177 + ,15.6 + ,3.5 + ,26.8590290272655 + ,11.6240545816589 + ,21.8017308111392 + ,10 + ,1.25 + ,24.5231261245390 + ,10.1198715180038 + ,3.80173081107135 + ,0 + ,3.0125 + ,22.1958135120851 + ,8.64677576944595 + ,6.30173081107452 + ,0.5 + ,0.5 + ,20.2774821608766 + ,7.20503203420426 + ,3.5517308110735 + ,0 + ,0 + ,18.2997339447889 + ,5.79327468630679 + ,-11.6982691889900 + ,0 + ,0.875 + ,18.2997339447889 + ,4.41123208898345 + ,-7.69826918897657 + ,0.75 + ,3.125 + ,15.0520550457364 + ,3.05873288600178 + ,31.5517308111828 + ,0.75 + ,10 + ,13.9578298121826 + ,1.73578547934020 + ,-2.44826918899326 + ,0 + ,0 + ,13.5917030950519 + ,0.443415523420818 + ,13.8017308111198 + ,0 + ,21 + ,13.5917030950519 + ,-0.822338412692918 + ,-12.9457691890859 + ,2.5 + ,0 + ,13.0645135486197 + ,-2.05516940597415 + ,-10.9482691889868 + ,0 + ,0.4125 + ,13.0645135486197 + ,-3.26375371120254 + ,-4.69826918893965 + ,0) + ,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) Warning message: package 'rpart' was built under R version 2.8.1 and help may not work correctly > 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.5689519 0 1.0000000 1.067115 0.3120007 2 0.0100000 1 0.4310481 0.897432 0.2268145 > postscript(file="/var/www/rcomp/tmp/1eqj51274872502.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/rcomp/tmp/2eqj51274872502.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 46.0000 27.500250 18.499750 2 40.5000 27.500250 12.999750 3 22.5000 27.500250 -5.000250 4 25.0000 27.500250 -2.500250 5 22.2500 27.500250 -5.250250 6 7.0000 4.628333 2.371667 7 11.0000 4.628333 6.371667 8 50.2500 27.500250 22.749750 9 16.2500 27.500250 -11.250250 10 32.5000 27.500250 4.999750 11 5.7525 27.500250 -21.747750 12 7.7500 4.628333 3.121667 13 14.0000 27.500250 -13.500250 14 3.5000 4.628333 -1.128333 15 1.2500 4.628333 -3.378333 16 3.0125 4.628333 -1.615833 17 0.5000 4.628333 -4.128333 18 0.0000 4.628333 -4.628333 19 0.8750 4.628333 -3.753333 20 3.1250 4.628333 -1.503333 21 10.0000 4.628333 5.371667 22 0.0000 4.628333 -4.628333 23 21.0000 4.628333 16.371667 24 0.0000 4.628333 -4.628333 25 0.4125 4.628333 -4.215833 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/rcomp/tmp/3eqj51274872502.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/421cd1274872502.tab") > > try(system("convert tmp/1eqj51274872502.ps tmp/1eqj51274872502.png",intern=TRUE)) character(0) > try(system("convert tmp/2eqj51274872502.ps tmp/2eqj51274872502.png",intern=TRUE)) character(0) > try(system("convert tmp/3eqj51274872502.ps tmp/3eqj51274872502.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.040 0.680 1.405