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(594.25 + ,NA + ,699.997920798106 + ,593.655750472382 + ,722.75 + ,853.75 + ,594.25 + ,648.094696588097 + ,647.942799207107 + ,803.8 + ,766.5 + ,620.2 + ,749.03449041482 + ,744.673364998473 + ,876.8 + ,758.05 + ,634.83 + ,757.606916738666 + ,753.778778169047 + ,795.8 + ,756.85 + ,647.152 + ,757.824390989537 + ,755.529720168062 + ,917.3 + ,685.4 + ,658.1218 + ,757.346140130168 + ,749.907388162832 + ,733 + ,696.525 + ,660.84962 + ,722.033515408683 + ,723.33892944098 + ,850.5 + ,610.025 + ,664.417158 + ,709.513418705798 + ,704.516968331494 + ,678.8 + ,708.325 + ,658.9779422 + ,660.682485072171 + ,672.68501267814 + ,910.3 + ,619.1 + ,663.91264798 + ,684.066397640522 + ,680.23484230233 + ,730 + ,740.525 + ,659.431383182 + ,652.179572089688 + ,664.582874665436 + ,798.5 + ,730.5 + ,667.5407448638 + ,695.541299848975 + ,696.052601976007 + ,763.5 + ,489.75 + ,673.83667037742 + ,712.699738883378 + ,655.544241875514 + ,725.3 + ,766.525 + ,655.428003339678 + ,603.271485377383 + ,680.45740717072 + ,768.5 + ,780.125 + ,666.53770300571 + ,683.399621312151 + ,686.76864206964 + ,NA + ,804.975 + ,677.896432705139 + ,730.874398867526 + ,724.218996193073 + ,883 + ,529.25 + ,690.604289434625 + ,767.244476743744 + ,753.807367605178 + ,2289 + ,743.75 + ,674.468860491163 + ,650.431960772313 + ,635.506534904222 + ,775.5 + ,771.15 + ,681.396974442046 + ,696.234347016108 + ,686.076054174035 + ,684 + ,830.5 + ,690.372276997842 + ,733.004468871136 + ,691.728771063574 + ,685.5 + ,600 + ,704.385049298058 + ,780.857252864086 + ,782.989097499081 + ,843.8 + ,856.1 + ,693.946544368252 + ,692.088845177539 + ,676.610491913217 + ,855.3 + ,702.75 + ,710.161889931427 + ,772.588846288154 + ,794.826530239656 + ,839.3 + ,533.775 + ,709.420700938284 + ,738.31052412253 + ,749.345042400237 + ,708.8 + ,311.25 + ,691.856130844456 + ,637.920340681366 + ,576.525816892812 + ,546.3 + ,590 + ,653.79551776001 + ,477.583912246844 + ,540.525418350128 + ,644.75 + ,738 + ,647.415965984009 + ,532.760008036036 + ,581.845867357353 + ,805.5 + ,797.05 + ,656.474369385608 + ,633.495958581145 + ,660.565830930517 + ,772.5 + ,531.3 + ,670.531932447048 + ,713.771599162515 + ,600.288183112223 + ,756.6 + ,820 + ,656.608739202343 + ,624.210837572008 + ,674.707766590396 + ,800 + ,533.25 + ,672.947865282109 + ,720.30812930781 + ,745.636209048273 + ,780 + ,633.25 + ,658.978078753898 + ,628.496205702196 + ,694.548612624309 + ,621 + ,634.275 + ,656.405270878508 + ,630.829464355652 + ,554.468246486220 + ,643 + ,747.3 + ,654.192243790657 + ,632.520603129091 + ,712.272349523156 + ,668 + ,220.375 + ,663.503019411592 + ,688.856658963743 + ,640.90531318958 + ,635 + ,195.75 + ,619.190217470432 + ,458.916359429194 + ,391.230261941304 + ,306 + ,123.25 + ,576.846195723389 + ,329.748972922716 + ,245.312950038669 + ,274 + ,161.75 + ,531.486576151050 + ,228.395088989555 + ,281.611155066594 + ,305 + ,126.75 + ,494.512918535945 + ,195.684327678938 + ,286.879143124124 + ,339 + ,285.1 + ,457.736626682351 + ,161.849961591333 + ,239.469364780464 + ,220 + ,461.5 + ,440.472964014116 + ,222.343580086496 + ,194.082609246879 + ,340 + ,463.625 + ,442.575667612704 + ,339.726401331110 + ,389.383530982008 + ,1651 + ,325.875 + ,444.680600851434 + ,400.538346355494 + ,290.573224811809 + ,448 + ,177 + ,432.800040766290 + ,363.892061710736 + ,332.245829765658 + ,243 + ,223 + ,407.220036689661 + ,272.161647449797 + ,316.860055686865 + ,211 + ,168.45 + ,388.798033020695 + ,248.032113830012 + ,272.113450680930 + ,190 + ,251.75 + ,366.763229718626 + ,208.971598350831 + ,33.9050653080551 + ,157 + ,131.5 + ,355.261906746763 + ,229.968105472534 + ,146.848855028763 + ,189 + ,110.375 + ,332.885716072087 + ,181.63796226685 + ,136.067930419851 + ,135 + ,164.125 + ,310.634644464878 + ,146.660655096268 + ,103.076560780786 + ,164) + ,dim=c(5 + ,50) + ,dimnames=list(c('Actuals' + ,'Croston' + ,'ETS' + ,'Arima' + ,'Kaneka') + ,1:50)) > y <- array(NA,dim=c(5,50),dimnames=list(c('Actuals','Croston','ETS','Arima','Kaneka'),1:50)) > 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.76661025 0 1.0000000 1.0251715 0.1286693 2 0.02004554 1 0.2333897 0.4141183 0.1026405 3 0.01000000 3 0.1932987 0.4118412 0.0850544 > postscript(file="/var/www/html/rcomp/tmp/11l0i1274871473.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/2cuh31274871473.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 594.250 613.6000 -19.3500000 2 853.750 765.8708 87.8791667 3 766.500 664.1143 102.3857143 4 758.050 664.1143 93.9357143 5 756.850 664.1143 92.7357143 6 685.400 664.1143 21.2857143 7 696.525 664.1143 32.4107143 8 610.025 664.1143 -54.0892857 9 708.325 765.8708 -57.5458333 10 619.100 765.8708 -146.7708333 11 740.525 765.8708 -25.3458333 12 730.500 765.8708 -35.3708333 13 489.750 664.1143 -174.3642857 14 766.525 765.8708 0.6541667 15 780.125 765.8708 14.2541667 16 804.975 664.1143 140.8607143 17 529.250 664.1143 -134.8642857 18 743.750 613.6000 130.1500000 19 771.150 765.8708 5.2791667 20 830.500 664.1143 166.3857143 21 600.000 664.1143 -64.1142857 22 856.100 765.8708 90.2291667 23 702.750 664.1143 38.6357143 24 533.775 664.1143 -130.3392857 25 311.250 239.4735 71.7764706 26 590.000 613.6000 -23.6000000 27 738.000 613.6000 124.4000000 28 797.050 765.8708 31.1791667 29 531.300 613.6000 -82.3000000 30 820.000 765.8708 54.1291667 31 533.250 664.1143 -130.8642857 32 633.250 239.4735 393.7764706 33 634.275 613.6000 20.6750000 34 747.300 765.8708 -18.5708333 35 220.375 239.4735 -19.0985294 36 195.750 239.4735 -43.7235294 37 123.250 239.4735 -116.2235294 38 161.750 239.4735 -77.7235294 39 126.750 239.4735 -112.7235294 40 285.100 239.4735 45.6264706 41 461.500 239.4735 222.0264706 42 463.625 613.6000 -149.9750000 43 325.875 239.4735 86.4014706 44 177.000 239.4735 -62.4735294 45 223.000 239.4735 -16.4735294 46 168.450 239.4735 -71.0235294 47 251.750 239.4735 12.2764706 48 131.500 239.4735 -107.9735294 49 110.375 239.4735 -129.0985294 50 164.125 239.4735 -75.3485294 > myr <- residuals(m) > myp <- predict(m) > postscript(file="/var/www/html/rcomp/tmp/3cuh31274871473.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/484xt1274871473.tab") > > try(system("convert tmp/11l0i1274871473.ps tmp/11l0i1274871473.png",intern=TRUE)) character(0) > try(system("convert tmp/2cuh31274871473.ps tmp/2cuh31274871473.png",intern=TRUE)) character(0) > try(system("convert tmp/3cuh31274871473.ps tmp/3cuh31274871473.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.073 0.474 1.252