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(250.75 + ,NA + ,250.850152884021 + ,251.27004 + ,336.675 + ,314.5125 + ,250.75 + ,314.369477869004 + ,251.27004 + ,352.25 + ,449.3885 + ,257.12625 + ,448.784077969190 + ,251.27004 + ,392.025 + ,305.7 + ,276.352475 + ,305.685475232751 + ,251.27004 + ,435.275 + ,162.375 + ,279.2872275 + ,162.837106781298 + ,251.27004 + ,148.75 + ,352.025 + ,267.59600475 + ,351.787128034443 + ,251.27004 + ,388.875 + ,379.125 + ,276.038904275 + ,378.830448248692 + ,251.27004 + ,359.65 + ,327.125 + ,286.3475138475 + ,327.047794937283 + ,251.27004 + ,139.9 + ,423.6625 + ,290.42526246275 + ,423.280779333534 + ,251.27004 + ,440.5 + ,152.25 + ,303.748986216475 + ,152.776138996337 + ,251.27004 + ,220.875 + ,183.8125 + ,288.599087594828 + ,184.199992170458 + ,251.27004 + ,234.5 + ,153.8875 + ,278.120428835345 + ,154.348519548268 + ,251.27004 + ,161.9 + ,245.625 + ,265.697135951810 + ,245.711928676033 + ,251.27004 + ,258.25 + ,108.9 + ,263.689922356629 + ,109.490720296634 + ,251.27004 + ,222.5 + ,291.625 + ,248.210930120966 + ,291.441298968799 + ,251.27004 + ,279.35 + ,284.875 + ,252.552337108870 + ,284.744734135499 + ,251.27004 + ,415.8 + ,192.25 + ,255.784603397983 + ,192.499491581822 + ,251.27004 + ,168.9 + ,45.2625 + ,249.431143058184 + ,46.0800968454689 + ,251.27004 + ,67.125 + ,205.375 + ,229.014278752366 + ,205.452810842879 + ,251.27004 + ,170.5 + ,301.25 + ,226.650350877129 + ,300.883977236240 + ,251.27004 + ,176.25 + ,165.375 + ,234.110315789416 + ,165.658013776605 + ,251.27004 + ,175.25 + ,281.6375 + ,227.236784210475 + ,281.371526167842 + ,251.27004 + ,262.625 + ,140.5875 + ,232.676855789427 + ,140.975949750512 + ,251.27004 + ,263.3 + ,331.75 + ,223.467920210485 + ,331.23300462894 + ,251.27004 + ,212.275 + ,232.625 + ,234.296128189436 + ,232.626038312863 + ,251.27004 + ,171.775) + ,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.7028891 0 1.0000000 1.0727678 0.2581946 2 0.0100000 1 0.2971109 0.4006243 0.1070606 > postscript(file="/var/www/html/freestat/rcomp/tmp/1zagu1274872684.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/2s1fx1274872684.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 250.7500 324.2179 -73.467929 2 314.5125 324.2179 -9.705429 3 449.3885 324.2179 125.170571 4 305.7000 324.2179 -18.517929 5 162.3750 158.4273 3.947727 6 352.0250 324.2179 27.807071 7 379.1250 324.2179 54.907071 8 327.1250 324.2179 2.907071 9 423.6625 324.2179 99.444571 10 152.2500 158.4273 -6.177273 11 183.8125 158.4273 25.385227 12 153.8875 158.4273 -4.539773 13 245.6250 324.2179 -78.592929 14 108.9000 158.4273 -49.527273 15 291.6250 324.2179 -32.592929 16 284.8750 324.2179 -39.342929 17 192.2500 158.4273 33.822727 18 45.2625 158.4273 -113.164773 19 205.3750 158.4273 46.947727 20 301.2500 324.2179 -22.967929 21 165.3750 158.4273 6.947727 22 281.6375 324.2179 -42.580429 23 140.5875 158.4273 -17.839773 24 331.7500 324.2179 7.532071 25 232.6250 158.4273 74.197727 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/freestat/rcomp/tmp/3s1fx1274872684.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/46bvo1274872684.tab") > > try(system("convert tmp/1zagu1274872684.ps tmp/1zagu1274872684.png",intern=TRUE)) character(0) > try(system("convert tmp/2s1fx1274872684.ps tmp/2s1fx1274872684.png",intern=TRUE)) character(0) > try(system("convert tmp/3s1fx1274872684.ps tmp/3s1fx1274872684.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.392 0.719 1.588