calif <- read.table("http://www.stat.cmu.edu/~cshalizi/350/hw/06/cadata.dat", header = TRUE)
### California Housing Dataset
price.deciles <- quantile(calif$MedianHouseValue, 0:10/10)
cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE)
plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20, xlab = "Longitude", ylab = "Latitude")
## gusca podrucja s visim cijenama su gradovi (LA, San Francisco...)

require(tree)
treefit <- tree(log(MedianHouseValue) ~ Longitude + Latitude, data = calif)

plot(treefit)
text(treefit, cex = 0.75) 
## u listovima pisu prosjecne (log) vrijednosti u odgovarajucem elementu particije
## za novu kombinaciju (Long, Lat) jednostavno se spustamo po stablu sve dok ne dodemo do lista
# te je vrijednost u tom listu upravo predikcija koju daje nas model
## velicina bridova odgovara smanjenju troska C(T) (vidi (7.4))


# buduci da imamo samo dva prediktora, mozemo prikazati rezultirajucu particiju prostora prediktora
cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE)
plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20,
     xlab = "Longitude", ylab = "Latitude")
partition.tree(treefit, ordvars = c("Longitude", "Latitude"), add = TRUE)


summary(treefit)
### Residual mean deviance (0.1662) je skoro trening MSE (C(T) u (7.3)) -- RSS se dijeli
# s (n-broj listova) umjesto n 

?tree ## kriterij zaustavljanja mozemo mijenjati koristeci parametre minsize, mincut i mindev
?tree.control

## smanjimo mindev -> dakle, napravit cemo neke podjele koje prvi model ne bi napravio
treefit2 <- tree(log(MedianHouseValue) ~ Longitude + Latitude, data = calif, mindev = 0.001)
plot(treefit2)
text(treefit2, cex = 0.75)

summary(treefit2) ## 68 listova

cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE)
plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20,
     xlab = "Longitude", ylab = "Latitude")
partition.tree(treefit2, ordvars = c("Longitude", "Latitude"), add = TRUE)
## puno finija podjela u gradovima

## ukljucimo i ostale prediktore
treefit3 <- tree(log(MedianHouseValue) ~ ., data = calif)

plot(treefit3)
text(treefit3, cex = 0.75)
## iako dopustamo samo binarne podjele (x_j < ili > od s), slozenije podjele mozemo dobiti
# koristeci nekoliko binarnih -> vidi MedianIncome (cini se da je to najznacajniji prediktor)
## stabla zapravo traze (moguce) komplicirane interakcije medu varijablama -> npr. Lat i Long se koriste
# samo ako je medianIncome<3.54, ali drugacije ovisno o tome je li < ili > od 2.51 itd. (npr. kod lin. regresije ili varijanti
# te interekcije ce se tesko prepoznati ako jer moramo ukljuciti puno varijabli)

summary(treefit3) ## 15 listova, te koristimo samo 4 od 8 prediktora
treefit3 ## detalji oko izbora podjela (deviance je suma kvadrata reziduala samo u tom cvoru)

## ne mozemo crtati dobivenu particiju, ali mozemo prikazati dobivene procjene kao sto smo prikazali stvarne cijene
cut.predictions <- cut(predict(treefit3), log(price.deciles), include.lowest = TRUE)
plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.predictions],
     pch = 20, xlab = "Longitude", ylab = "Latitude")

cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE)
plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20,
     xlab = "Longitude", ylab = "Latitude")

## obrezivanje i CV
?cv.tree
?prune.tree  ## daje niz podstabala iz (7.9)

pr.tree=prune.tree(treefit3)
pr.tree  
## size su njihove velicine, dev je RSS, a k je zapravo alfa

pr.tree=prune.tree(treefit3, best=8) ## probaj 8
plot(pr.tree)
text(pr.tree)
summary(pr.tree)

cv.calif=cv.tree(treefit3)
names(cv.calif)
plot(cv.calif, type="b")
## u ovom slucaju, najbolji izgleda puni model 

### pogledajmo kakva je cv procjena testne greske u odnosu na stvarnu

set.seed(1021)
n=20640
m=10000 ## velicina skupa za trening
train=sample(1:n, m)
treefit4 <- tree(log(MedianHouseValue) ~ ., data = calif, subset=train)
summary(treefit4)
cv.calif=cv.tree(treefit4)
plot(cv.calif, type="b")
which.min(cv.calif$dev)

y.test=log(calif$MedianHouseValue[-train])
test.er=numeric(14)
for (i in 2:14) {
  tr=prune.tree(treefit4, best = i)
  yhat=predict(tr, newdata=calif[-train, ])
  test.er[i]=mean((yhat-y.test)^2)
}

plot(cv.calif$size, cv.calif$dev/m, type="b")  ## dev/m=RSS/m=MSE   
points(2:14, test.er[-1], type="b", col="blue")

## zbog ogromog broja podataka, tesko cemo imati overfitting

## testna greska punog modela (sa 14 cvorova)
test.er[14]
## [1] 0.1404545 (gledamo logaritme cijena pa je stoga greska manjeg reda)

######## Random forrest

library(randomForest)

rf.cal<- randomForest(log(MedianHouseValue) ~ ., data = calif,
                        subset = train, importance = TRUE, ntree=1000)

## treba malo duze vremena da se izvrti

rf.cal  # m=2

plot(1:1000, rf.cal$mse, type="l")
yhat=predict(rf.cal, newdata=calif[-train, ])
mean((yhat-y.test)^2)
## [1] 0.062861  -- puno bolje nego jedno stablo

varImpPlot(rf.cal, type=2) ## slicno kao za Boston podatke

## jedno stablo:
summary(treefit4)
plot(treefit4)

text(treefit4, cex = 0.75)