# Lab: Classification Methods


## The Stock Market Data

###
library(ISLR2)
?Smarket

head(Smarket)
names(Smarket)
dim(Smarket)
summary(Smarket)

###
pairs(Smarket)
# Direction == "up" akko Today>0
cor(Smarket[, -9])
# velika korelacija izmedu Volume i Year; mala korelacija Today i Lag_k

boxplot(Volume~Year, data=Smarket)


attach(Smarket)
###
train <- (Year < 2005)
Smarket.2005 <- Smarket[!train, ]
dim(Smarket.2005)
# 252 povrata u 2005
n=1250-252 ## velicina skupa za trening
Direction.2005 <- Direction[!train]

############################ logisticka regresija ###################################

# Neka je Y=Direction te pretpostavljamo model 
# P(Y="Up" | Lag1,...., Lag5, Volume) = h(b0 + b1*Lag1 + ... + b5*Lag5 + b6*Volume), za
# h(z) = exp(z)/(1+exp(z))

glm.fits <- glm(
  Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
  data = Smarket, family = binomial, subset = train
)


summary(glm.fits)
## procijenjeni parametri (b0_hat,..., b6_hat) su dakle
coef(glm.fits)

# p-vrijednost u zadnjem stupcu odnosi se na testiranje H_0: bj=0 pri cemu je testna
# statistika slicna kao kod lin. regresije (tzv. Waldova statistika), te rezultate slicno intepretiramo
# -> niti jedna se zasebno ne cini znacajna (najmanja p-vrijednost za Lag1)

# glm.fits0 <- glm(
#   Direction ~ 1,
#   data = Smarket, family = binomial, subset = train
# )
# 
# anova(glm.fits0, glm.fits, test="chisq")
# ### velika p-vrijednost --  ne mozemo odbaciti nultu hipotezu da je model samo sa slobodnim clanom
# ### (dakle, bez kovarijata) dovoljan.


## procjena testne greske
glm.probs <- predict(glm.fits, Smarket.2005,
                     type = "response")

glm.pred <- rep("Down", 252)
glm.pred[glm.probs > .5] <- "Up"

table(glm.pred, Direction.2005) # kontingencijska tablica za vektor (glm.pred, Direction.2005)
#          Direction.2005
# glm.pred Down Up
# Down       77 97
# Up         34 44
### procjena testne greske P(Y_hat != Y) (dakle, uz 0-1 funkciju gubitka) je
mean(glm.pred != Direction.2005)
##  0.5198413 = (34+97)/252 --> losije od slucajnog pogadanja!

## probajmo:
glm.fits <- glm(Direction ~ Lag1 + Lag2, data = Smarket,
                family = binomial, subset = train)
glm.probs <- predict(glm.fits, Smarket.2005,
                     type = "response")
glm.pred <- rep("Down", 252)
glm.pred[glm.probs > .5] <- "Up"
table(glm.pred, Direction.2005)
#          Direction.2005
# glm.pred Down  Up
# Down       35  35
# Up         76 106

mean(glm.pred != Direction.2005)
# [1] 0.4404762 --> izgleda puno bolje, ali
# kada bi stavili f_hat(x):= "up" za sve x, dobili bi istz testnu gresku (to je vjerojatno slucajnost)
mean(Direction.2005 != "Up")
# [1] 0.4404762

# Ipak, uocimo da je u nasem modelu procjena za uvjetnu gresku P(Y != Y_hat | Y_hat="Up")
76/(76+106)
# [1] 0.4175824

## Primjer predikcije

h=function(z) exp(z)/(1+exp(z))

# nas model je sada 
# P(Y="Up" | Lag1, Lag2) = h(b0 + b1*Lag1 + b2*Lag2)

b=coef(glm.fits)
# procjena za P(Y="Up" | Lag1=1.2, Lag2=1.1) je
h(b[1] + b[2]*1.2 + b[3]*1.1)
# 0.4791462 

# ili preko predict:
predict(glm.fits,
        newdata =
          data.frame(Lag1 = 1.2,  Lag2 = 1.1),
        type = "response"
)
# 0.4791462 


######################################### LDA ###############################

library(MASS)
lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket,
               subset = train)
# sintaksa je ista kao kod glm() funkcije, samo bez opcije "family"

lda.fit
# dakle, ako je 0="Down", 1="Up", 
# 1. pi0_hat=0.491984, pi1_hat=0.508016
# 2. mi0_hat=(0.04279022, 0.03389409) approx E(X | Y=0), mi1_hat=(-0.03954635 -0.03132544)
# 3. Coefficients of linear discriminants -- o tome cemo malo kasnije

lda.pred <- predict(lda.fit, Smarket.2005)
names(lda.pred)

# lda.pred$posterior daje P(Y= "Down" | X=x) i P(Y= "Up" | X=x)

lda.class <- lda.pred$class
table(lda.class, Direction.2005)

mean(lda.class != Direction.2005)
# [1] 0.4404762 -- isto kao kod log. regresije jer daju potpuno iste predikcije na testnom skupu:

table(lda.class, glm.pred)



# ### Coefficients of linear discriminants:
# w=as.vector(lda.fit$scaling) ## <-- vektor koeficijenata (ili smjer) prve disktriminacijske koordinate (jedina jer K=1)
# pi=as.vector(lda.fit$prior)
# mi_0=as.vector(lda.fit$means[1,])
# mi_1=as.vector(lda.fit$means[2,])
# 
# lda.pred <- predict(lda.fit, Smarket.2005)$class 
# # daje procijenjene klase na testnom skupu
# 
# ### pokazimo da je 
# ### f_hat(x) = argmax_{k=0,1} delta_k(x) = argmin_{k=0,1} [ ||w*x - w*mi_k_hat ||/2 - log(pi_k_hat) ]
# delta_0=function(x){abs(w%*%x - w%*%mi_0)^2/2 - log(pi[1])}
# delta_1=function(x){abs(w%*%x - w%*%mi_1)^2/2 - log(pi[2])}
# f_hat=function(x){
#   if (delta_0(x)<delta_1(x)) {
#     return("Down") 
#   }
#   else {return("Up")}
# }
# 
# data=as.matrix(Smarket.2005[,c(2,3)])
# lda.pred_2 = apply(data, 1, f_hat)
# table(lda.pred, lda.pred_2)


########################## QDA ########################################
qda.fit <- qda(Direction ~ Lag1 + Lag2, data = Smarket,
               subset = train)
qda.fit

qda.class <- predict(qda.fit, Smarket.2005)$class
table(qda.class, Direction.2005)

mean(qda.class != Direction.2005)
# [1] 0.4007937 -- puno bolje nego log. regr. i LDA

######################## Naive Bayes #####################
library(e1071)
nb.fit <- naiveBayes(Direction ~ Lag1 + Lag2, data = Smarket,
                     subset = train)
### pretpostavlja se da su uvjetno na Y=k, kovarijate nezavisne i normalno distribuirane
nb.fit
### npr. uvjetno na Y="Down", procijenjeno ocekivanje i varijanca za Lag1 su 0.04279022 i 1.227446.

nb.class <- predict(nb.fit, Smarket.2005)
table(nb.class, Direction.2005)

mean(nb.class != Direction.2005)
## [1] 0.4087302 -- samo malo losije od QDA, uocimo:
table(qda.class, nb.class)

## izgleda da LDA i log. regresija u ovom slucaju imaju preveliku pristranost jer granica odluke nije priblizno linearna

## procjena vjerojatnosti P(Y=k | X=x) na testnom skupu
nb.preds <- predict(nb.fit, Smarket.2005, type = "raw")
nb.preds[1:5, ]


#######################  kNN ################################
 
#sintaksa je malo drugacija

library(class)
train.X <- cbind(Lag1, Lag2)[train, ]
test.X <- cbind(Lag1, Lag2)[!train, ]
train.Direction <- Direction[train]

set.seed(1) ## kako bi mogli reproducirati rezultate jer
             # u slucaju da ima vise podatkaka koji su k-ti susjed, slucajno se odabire jedan
knn.pred <- knn(train.X, test.X, train.Direction, k = 1)
table(knn.pred, Direction.2005)

mean(knn.pred != Direction.2005)
# [1] 0.5 -- k=1 daje preveliku varijancu

set.seed(5)
knn.pred=numeric(300)
for (i in 1:300) {
  pred <- knn(train.X, test.X, train.Direction, k = i)
  knn.pred[i]=mean(pred != Direction.2005)
}
c(which.min(knn.pred), min(knn.pred))

### [1] 242.0000000   0.3888889 -- izgleda super! U cemu je problem?

### -> Ovo nije dobra procjena za testnu gresku!
### Opcenito, skup za odabir parametra kompleksnosti i procjenu greske ne smije biti isti! 
### Koristimo dodatni skup za tzv. validaciju.