# początek podobny jak poprzednio
library(mlbench)
data(HouseVotes84)
library(rpart)

      # podzielmy dane na trenujące i testowe losowo w stosunku ok.  2:1
rhv <- runif(nrow(HouseVotes84))
hv.train <- HouseVotes84[rhv>=0.33,]
hv.test <- HouseVotes84[rhv<0.33,]

      # pakiet rpart zawiera implementację indukcji drzew decyzyjnych
library(rpart)

      # drzewo dla domyślnych parametrów
hv.tree0 <- rpart(Class ~ ., hv.train)
      # drzewo z maksymalnie rozluźnionym kryterium stopu
hv.tree1 <- rpart(Class ~ ., hv.train, minsplit=2, cp=0)

      # wykaz możliwości przycięcia dla różnych wartości cp
hv.tree1$cptable
      # przykładowe drzewa uzyskiwane w wyniku przycięcia
hv.tree2 <- prune(hv.tree1, cp=0.001)
hv.tree3 <- prune(hv.tree1, cp=0.005)
hv.tree4 <- prune(hv.tree1, cp=0.01)

      # testowanie
err <- function(y.true, y.pred) { sum(y.pred!=y.true)/length(y.true) }
err(hv.test$Class, predict(hv.tree1, hv.test, type="c"))
err(hv.test$Class, predict(hv.tree2, hv.test, type="c"))
err(hv.test$Class, predict(hv.tree3, hv.test, type="c"))
err(hv.test$Class, predict(hv.tree4, hv.test, type="c"))

      # sprawdzamy dostępność pakietu e1071 i instalujemy w razie potrzeby
if (! "e1071" %in% row.names(installed.packages()))
  install.packages("e1071")
     # ładujemy pakiet w celu użycia funkcji tune
library(e1071)
  
      # poniższa linijka potrzebna do ominięcia błędu w funkcji tune
attr(hv.train, "na.action") <- na.rpart
      # potrzebne jest też opakowanie na funkcję predict
      # niewymagające przekazywania argumentu type
predc <- function(...) { predict(..., type="c") }

      # 10-krotna walidacja krzyżowa na zbiorze trenującym
      # przy domyślnych parametrach
hv.cv10 <- tune(rpart, Class ~ ., data=hv.train, na.action=na.rpart,
                predict.func=predc,
                tunecontrol=tune.control(sampling="cross", cross=10))
      # sprawdźmy błąd uzyskany dla walidacji krzyżowej
hv.cv10
      # końcowy model zbudowanego na całości przekazanych danych
hv.cv10$best.model
      # to powinno być to samo co wcześniej przy bezpośrednim wywołaniu rpart
hv.tree0
      # zobaczmy jak błąd na zbiorze testowym ma się do błędu z walidacji krzyżowej
err(hv.test$Class, predict(hv.cv10$best.model, hv.test, type="c"))

      # strojenie parametrów
      # badane wartości dla minsplit, cp, usesurrogate
par.minsplit <- c(2, 5, 10, 20, 50)
par.cp <- c(0, 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5)
par.usesurrogate <- c(0, 1, 2)
      # pominiemy argument tunecontrol, gdyż i tak
      # 10-krotna walidacja krzyżowa używana jest domyślnie
      # uwaga -- to może trochę potrwać
hv.tune <- tune(rpart, Class ~ ., data=hv.train, na.action=na.rpart,
                ranges=list(minsplit=par.minsplit, cp=par.cp,
                            usesurrogate=par.usesurrogate),
                predict.func=predc)
      # wynik strojenia
hv.tune
      # najlepszy zestaw parametrów
hv.tune$best.parameters
      # model dla najlepszego zestawu parametrów zbudowany
      # na całym przekazanym zbiorze danych
hv.tune$best.model
      # czy i ewentualnie czym to się różni od drzewa
      # dla parametrów domyślnych?
summary(hv.tree0)
summary(hv.tree.best)
      # w szczególności, zobaczmy, które przykłady trafiają
      # do innych liści
hv.train[hv.tree0$where!=hv.tree.best$where,]
      # wyniki dla wszystkich badanych zestawów parametrów
      # uporządkowane wg błędu
hv.tune$performances[order(hv.tune$performances$error),]
      # ponownie zbudujmy model dla najlepszego zestawu
hv.tree.best <- rpart(Class ~ ., hv.train,
                      control=rpart.control(hv.tune$best.parameters))
      # czy wyszło to samo drzewo, które jest zachowane jako hv.tune$best.model?
hv.tree.best
      # sprawdźmy błąd na zbiorze testowym
err(hv.test$Class, predict(hv.tree.best, hv.test, type="c"))