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(719.625 + ,NA + ,712.710818344606 + ,718.905375359812 + ,619.5 + ,689.5125 + ,719.625 + ,718.479353436485 + ,719.62500000013 + ,659.125 + ,640.5 + ,716.61375 + ,694.312167950503 + ,689.5125 + ,675.65 + ,541.3875 + ,709.002375 + ,649.416416773587 + ,640.5 + ,555.5 + ,521.3875 + ,692.2408875 + ,559.287370647071 + ,541.3875 + ,582.75 + ,613.325 + ,675.15554875 + ,527.667325831681 + ,521.3875 + ,466.75 + ,683 + ,668.972493875 + ,599.131936389492 + ,613.325 + ,562.8 + ,728.6375 + ,670.3752444875 + ,669.103466230299 + ,683 + ,289.65 + ,452.5375 + ,676.20147003875 + ,718.772997648195 + ,728.6375 + ,592 + ,380.425 + ,653.835073034875 + ,496.65143830366 + ,452.5375 + ,461.25 + ,377.125 + ,626.494065731387 + ,399.683160440186 + ,380.425 + ,429.15 + ,316.75 + ,601.557159158249 + ,380.862778420582 + ,377.125 + ,345.9 + ,387.2625 + ,573.076443242424 + ,327.373178264001 + ,316.75 + ,421 + ,409.75 + ,554.495048918181 + ,377.339128145203 + ,387.2625 + ,375.9 + ,497.875 + ,540.020544026363 + ,404.379674802946 + ,409.75 + ,438.4 + ,616.4 + ,535.805989623727 + ,482.383275403329 + ,497.875 + ,562.4 + ,715.5125 + ,543.865390661354 + ,594.194074897074 + ,616.4 + ,700.65 + ,454.925 + ,561.030101595219 + ,695.410621898232 + ,715.5125 + ,547.525 + ,464.25 + ,550.419591435697 + ,494.772308045129 + ,454.925 + ,484.5 + ,247.675 + ,541.802632292127 + ,469.307399279518 + ,464.25 + ,357.65 + ,292.5125 + ,512.389869062915 + ,284.398419958176 + ,247.675 + ,298 + ,416.525 + ,490.402132156623 + ,291.168036059174 + ,292.5125 + ,246.875 + ,481.6625 + ,483.014418940961 + ,395.75395548464 + ,416.525 + ,609.625 + ,219.7625 + ,482.879227046865 + ,467.427868382560 + ,481.6625 + ,359.9 + ,402.625 + ,456.567554342178 + ,260.799457420375 + ,219.7625 + ,402.875) + ,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.6393643 0 1.0000000 1.0713938 0.2140061 2 0.0100000 1 0.3606357 0.7668597 0.2994565 > postscript(file="/var/www/html/freestat/rcomp/tmp/16i891274873533.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/26i891274873533.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 719.6250 614.7292 104.8958333 2 689.5125 614.7292 74.7833333 3 640.5000 614.7292 25.7708333 4 541.3875 614.7292 -73.3416667 5 521.3875 614.7292 -93.3416667 6 613.3250 614.7292 -1.4041667 7 683.0000 614.7292 68.2708333 8 728.6375 614.7292 113.9083333 9 452.5375 614.7292 -162.1916667 10 380.4250 376.4769 3.9480769 11 377.1250 376.4769 0.6480769 12 316.7500 376.4769 -59.7269231 13 387.2625 376.4769 10.7855769 14 409.7500 376.4769 33.2730769 15 497.8750 376.4769 121.3980769 16 616.4000 614.7292 1.6708333 17 715.5125 614.7292 100.7833333 18 454.9250 614.7292 -159.8041667 19 464.2500 376.4769 87.7730769 20 247.6750 376.4769 -128.8019231 21 292.5125 376.4769 -83.9644231 22 416.5250 376.4769 40.0480769 23 481.6625 376.4769 105.1855769 24 219.7625 376.4769 -156.7144231 25 402.6250 376.4769 26.1480769 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/freestat/rcomp/tmp/3hr8u1274873533.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/4la6z1274873533.tab") > > try(system("convert tmp/16i891274873533.ps tmp/16i891274873533.png",intern=TRUE)) character(0) > try(system("convert tmp/26i891274873533.ps tmp/26i891274873533.png",intern=TRUE)) character(0) > try(system("convert tmp/3hr8u1274873533.ps tmp/3hr8u1274873533.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.363 0.667 1.534