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(724 + ,NA + ,726.178301126702 + ,723.250932184245 + ,763.275 + ,762.275 + ,724 + ,724.988015993228 + ,708.951659916233 + ,836.3 + ,721.125 + ,727.8275 + ,745.362675029711 + ,709.122505687033 + ,825.15 + ,653.275 + ,727.15725 + ,732.118526227605 + ,688.027750893318 + ,764.65 + ,663.7125 + ,719.769025 + ,689.036201204222 + ,650.13316612633 + ,820.15 + ,735.5125 + ,714.1633725 + ,675.198617089505 + ,630.378745649382 + ,781 + ,628.1375 + ,716.29828525 + ,708.15582272578 + ,643.244311473 + ,746.9 + ,792.55 + ,707.482206725 + ,664.431555782846 + ,612.847273246156 + ,441.5 + ,636.5 + ,715.9889860525 + ,734.439084864176 + ,652.317786039178 + ,1532.25 + ,800.825 + ,708.04008744725 + ,680.922408372177 + ,621.589491648231 + ,684.75 + ,728.05 + ,717.318578702525 + ,746.440564025738 + ,660.842815301877 + ,849.55 + ,618.2625 + ,718.391720832272 + ,736.39144148354 + ,659.888667071522 + ,774.05 + ,450.625 + ,708.378798749045 + ,671.84245826112 + ,619.88373918091 + ,595.525 + ,767.525 + ,682.60341887414 + ,550.963003806743 + ,534.087282375116 + ,789 + ,675.65 + ,691.095576986727 + ,669.298582567855 + ,592.774139057074 + ,778.45 + ,583.25 + ,689.551019288054 + ,672.769176074555 + ,597.440724777523 + ,700.5 + ,690.7875 + ,678.920917359249 + ,623.853375069715 + ,567.280936817024 + ,655.65 + ,208.0625 + ,680.107575623324 + ,660.4280675445 + ,586.525266208221 + ,470.65 + ,142.5 + ,632.903068060991 + ,413.24277101072 + ,425.668945310045 + ,289.375 + ,205.925 + ,583.862761254892 + ,265.301289743849 + ,299.002937771023 + ,279.275 + ,462.5625 + ,546.068985129403 + ,232.856411440937 + ,240.53933481731 + ,995.4 + ,251.4375 + ,537.718336616463 + ,358.374292751612 + ,295.130267176479 + ,345.5 + ,195.725 + ,509.090252954817 + ,299.941015059498 + ,254.38547660905 + ,200.125 + ,191.625 + ,477.753727659335 + ,242.994446986505 + ,208.270449663097 + ,172.95 + ,137.25 + ,449.140854893401 + ,214.924733247989 + ,177.229915639438 + ,149.3) + ,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.6955167 0 1.0000000 1.0915772 0.1909219 2 0.0100000 1 0.3044833 0.9383504 0.3513416 > postscript(file="/var/www/html/rcomp/tmp/1ksx01274874036.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/2ksx01274874036.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 724.0000 657.7847 66.215278 2 762.2750 657.7847 104.490278 3 721.1250 657.7847 63.340278 4 653.2750 657.7847 -4.509722 5 663.7125 657.7847 5.927778 6 735.5125 657.7847 77.727778 7 628.1375 657.7847 -29.647222 8 792.5500 657.7847 134.765278 9 636.5000 657.7847 -21.284722 10 800.8250 657.7847 143.040278 11 728.0500 657.7847 70.265278 12 618.2625 657.7847 -39.522222 13 450.6250 657.7847 -207.159722 14 767.5250 657.7847 109.740278 15 675.6500 657.7847 17.865278 16 583.2500 657.7847 -74.534722 17 690.7875 657.7847 33.002778 18 208.0625 657.7847 -449.722222 19 142.5000 226.7179 -84.217857 20 205.9250 226.7179 -20.792857 21 462.5625 226.7179 235.844643 22 251.4375 226.7179 24.719643 23 195.7250 226.7179 -30.992857 24 191.6250 226.7179 -35.092857 25 137.2500 226.7179 -89.467857 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/rcomp/tmp/3vke31274874036.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/4y2cr1274874036.tab") > > try(system("convert tmp/1ksx01274874036.ps tmp/1ksx01274874036.png",intern=TRUE)) character(0) > try(system("convert tmp/2ksx01274874036.ps tmp/2ksx01274874036.png",intern=TRUE)) character(0) > try(system("convert tmp/3vke31274874036.ps tmp/3vke31274874036.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.924 0.473 1.268