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(66 + ,NA + ,65.9293366525043 + ,65.9340000655022 + ,55 + ,71 + ,66 + ,70.8680320462113 + ,67.4510134015007 + ,77 + ,50 + ,66.5 + ,50.2163309228226 + ,65.1096772035754 + ,66.4 + ,66 + ,64.85 + ,65.9827164753333 + ,63.3900574001068 + ,66 + ,66 + ,64.965 + ,65.9993404862515 + ,65.7651067120978 + ,55 + ,44 + ,65.0685 + ,44.3436189251445 + ,56.7231840728576 + ,44 + ,76.75 + ,62.96165 + ,76.6231024714092 + ,61.0647320143803 + ,55 + ,66 + ,64.340485 + ,66.0460290529034 + ,64.1024412875046 + ,33 + ,65.75 + ,64.5064365 + ,65.8141920877573 + ,55.3499659970663 + ,77 + ,76.75 + ,64.63079285 + ,76.674093729902 + ,72.1767781103031 + ,66 + ,65 + ,65.842713565 + ,65.1023563693362 + ,68.3625842416714 + ,77 + ,76 + ,65.7584422085 + ,75.9657948158125 + ,67.7363072850427 + ,66 + ,88 + ,66.78259798765 + ,87.819777379446 + ,74.280295689305 + ,88 + ,75.5 + ,68.904338188885 + ,75.5018569541544 + ,72.3141534058821 + ,88 + ,97.5 + ,69.5639043699965 + ,97.2291009532818 + ,78.2382585173593 + ,77.75 + ,98 + ,72.3575139329969 + ,97.7403340016694 + ,87.6827804286996 + ,88 + ,88.0115 + ,74.9217625396972 + ,87.8957556168923 + ,84.8565432816706 + ,99 + ,55.25 + ,76.2307362857275 + ,55.5590661695737 + ,95.469375671624 + ,77 + ,88.25 + ,74.1326626571547 + ,88.1610477971363 + ,86.4534430561582 + ,88 + ,87 + ,75.5443963914393 + ,86.940601684571 + ,81.6480522907467 + ,120 + ,75.5 + ,76.6899567522953 + ,75.592583208209 + ,68.8361541965404 + ,88 + ,88 + ,76.5709610770658 + ,87.955890234155 + ,85.1449520700068 + ,77 + ,74.5 + ,77.7138649693592 + ,74.6272559913572 + ,84.9718409849175 + ,101.25 + ,63.5 + ,77.3924784724233 + ,63.7653950698418 + ,77.7240276566534 + ,77 + ,94 + ,76.003230625181 + ,93.9259056838178 + ,79.8596794457166 + ,88) + ,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.6676578 0 1.0000000 1.0940291 0.2612326 2 0.0100000 1 0.3323422 0.5318093 0.1255753 > postscript(file="/var/www/html/freestat/rcomp/tmp/14fqs1274872994.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/24fqs1274872994.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 66.0000 66.67647 -0.6764706 2 71.0000 66.67647 4.3235294 3 50.0000 66.67647 -16.6764706 4 66.0000 66.67647 -0.6764706 5 66.0000 66.67647 -0.6764706 6 44.0000 66.67647 -22.6764706 7 76.7500 66.67647 10.0735294 8 66.0000 66.67647 -0.6764706 9 65.7500 66.67647 -0.9264706 10 76.7500 66.67647 10.0735294 11 65.0000 66.67647 -1.6764706 12 76.0000 66.67647 9.3235294 13 88.0000 91.09519 -3.0951875 14 75.5000 66.67647 8.8235294 15 97.5000 91.09519 6.4048125 16 98.0000 91.09519 6.9048125 17 88.0115 91.09519 -3.0836875 18 55.2500 66.67647 -11.4264706 19 88.2500 91.09519 -2.8451875 20 87.0000 91.09519 -4.0951875 21 75.5000 66.67647 8.8235294 22 88.0000 91.09519 -3.0951875 23 74.5000 66.67647 7.8235294 24 63.5000 66.67647 -3.1764706 25 94.0000 91.09519 2.9048125 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/freestat/rcomp/tmp/3w7pv1274872994.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/4ip611274872994.tab") > > try(system("convert tmp/14fqs1274872994.ps tmp/14fqs1274872994.png",intern=TRUE)) character(0) > try(system("convert tmp/24fqs1274872994.ps tmp/24fqs1274872994.png",intern=TRUE)) character(0) > try(system("convert tmp/3w7pv1274872994.ps tmp/3w7pv1274872994.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.383 0.689 1.549