### preuzeto od Cosma Shalizi, kolegij na CMU ###
### California Housing Dataset
## podaci od kucanstvima za razne dijelove Kalifornije
## podjela teritorija na 20640 dijelova podjednake populacije
calif=read.table("cadata.dat", header=TRUE)
names(calif)
head(calif)

### crtamo MedianHouseValue (koristeci skalu sivih boja) 
##  u odnosu na lokaciju
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)
?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 (8.4))
##  ---> dulji brid predstavlja vece smanjenje!

### 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 (8.3)) -- RSS se dijeli
## sa (n-broj listova) umjesto n 

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

###  smanjimo mindev (donja ograda na smanjenje troska)-> 
##  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 interakcije ce se tesko prepoznati 
## jer moramo unaprijed odrediti koje interakcije dopustamo 

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 (8.9)

pr.tree=prune.tree(treefit3)
pr.tree  
## size su njihove velicine (tj. broj listova), 
## 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")
## slicno


### 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 forest
# 
# 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)