Classification

Задача 1

\[ p(X) = \frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}}\\ 1-p(x) = 1 - \frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}} = \frac{1+e^{\beta_0+\beta_1X}-e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}}=\frac{1}{1+e^{\beta_0+\beta_1X}}\\ \frac{p(X)}{1-p(x)}=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}} * \frac{1+e^{\beta_0+\beta_1X}}{1}=e^{\beta_0+\beta_1X} \]

Задача 2

Именителот можеме слободно да го игнорираме, бидејќи е ист за секоја класа. Веднаш ќе преминеме на запишување на остатокот од изразот во логаритамска форма. Изразите А, B, C, … имаат максимум за исто к (добиени се со отстранување на собироците во кои не се појавува к). Претпоставуваме дека работиме со природен логаритам т.е. log(e) = 1.

\[ A=log(\pi_k)-log(\sqrt{2\pi\sigma})-(\frac{(x-\mu_k)^2}{2\sigma^2})log(e)\\ B=log(\pi_k)-\frac{(x-\mu_k)^2}{2\sigma^2}=log(\pi_k)-\frac{x^2+2x\mu_k-\mu_k^2}{2\sigma^2}\\ C=log(\pi_k)-\frac{x\mu_k}{\sigma^2}+\frac{\mu_k^2}{2\sigma^2}=\delta_k(x) \]

Задача 3

Слично на задача 2, изведуваме делта к за x, но овој пат со различна сигма за секоја класа. Од тоа што x на квадрат фигурира во делта к за x, следи дека бајесовиот класификатор е квадратен.

\[ A=log(\pi_k)-log(\sqrt{2\pi\sigma_k})-(\frac{(x-\mu_k)^2}{2\sigma_k^2})log(e)\\ B=log(\pi_k)-log(\sqrt{2\pi\sigma_k})-\frac{(x-\mu_k)^2}{2\sigma_k^2}=log(\pi_k)-log(\sqrt{2\pi\sigma_k})-\frac{x^2+2x\mu_k-\mu_k^2}{2\sigma_k^2}\\ C=log(\pi_k)log(\sqrt{2\pi\sigma_k})-\frac{x^2}{2\sigma_k^2}-\frac{x\mu_k}{\sigma_k^2}+\frac{\mu_k^2}{2\sigma_k^2}=\delta_k(x) \]

Задача 4

(a)

\[ 2\int_{0}^{0.05}(x + 0.05)dx + \int_{0.05}^{0.95}0.1dx=2\int_{0}^{0.05}xdx + \int_{0}^{0.05}0.1dx + \int_{0.05}^{0.95}0.1dx=\\ 0.05^2+0.095=0.0975 = 9.75\% \]

(b)

Решението е двоен интеграл на делови (слично како под а), со тоа што внатрешниот интеграл е горниот резултат. За поедноставен запис, внатрешниот интеграл ќе го третираме како пресметан.

\[ 2\int_{0}^{0.05}(x + 0.05) * 0.0975dx + \int_{0.05}^{0.95}0.1 * 0.0975dx=2\int_{0}^{0.05}xdx + \int_{0}^{0.05}0.1dx + \int_{0.05}^{0.95}0.1dx=\\ 0.0975 * (0.05^2+0.095)=0.0975^2 = 0.00950625 = 0.950625\% \]

(c)

Шемата е проста. Со секое наголемување на димензијата, процентот се квадрира. За 100 добиваме:

\[ 0.0975^{2^{99}} \approx 0 \]

(d)

Очигледно е дека со наголемување на димензиите, двојно експоненцијално расте бројот на неискористени примероци за класификација. Многу брзо бројот на блиски премероци станува многу низок.

(e)

За p=1, 10% од должината на линијата е 0.1, па страната на хиперкоцката со димензија 1 е 0.1.

За p=2, 10% од плоштината на квадратот 0.1, па страната на хиперкоцката е sqrt(0.1) = 0.3162278

За p=100, 10% од мерата на хиперповршината е 0.1, па страната на хиперкоцката е а, каде што а^100=0.1.

Задача 5

(a)

LDA, и на тренирачкото и на тестирачкото множество.

(b)

QDA, и на тренирачкото и на тестирачкото множество.

(c)

Во општ случај, QDA треба да станува подобро со зголемување на бројот на примероци, бидејќи ќе се намалува грешката поради степенот на over-fitting на QDA, додека грешката поради степенот на under-fitting на LDA нема многу да се промени. Единствено кога бајесовата decision boundary е линеарна, одговорот ќе биде поинаков.

(d)

Не. QDA во најдобар случај ќе направи over-fitting на тренирачкото множество.

Задача 6

(a)

e = exp(-6 + 0.05 * 40 + 1 * 3.5)
probability = e / (1 + e)
print(probability)

[1] 0.3775407

(b)

desiredLogOdds = 1
currentLogOdds = probability / (1 - probability)
additionalHours = (desiredLogOdds - currentLogOdds) / 0.05
print(additionalHours)

[1] 7.869387

Задача 7

means = c(10, 0)
sds = c(6, 6)
pdfsTimesPriors = c(0, 0)
priors = c(0.8, 0.2)
for (k in 1:2) {
  pdfsTimesPriors[k] = priors[k] * pnorm(4, means[k], sds[k])
}
print(pdfsTimesPriors[1] / sum(pdfsTimesPriors))

[1] 0.4591621

Задача 8

Бидејќи користиме KNN со к=1, error rate на тренирачкото множество ќе ни биде 0 (на секоја точка, најблиската е самата таа и само таа ја земаме предвид). Бидејќи просекот е 18% error rate, error rate на тестирачкото множество е 36%, што е полошо од 30% на логистичка регресија. Бираме логистичка регресија.

Задача 9

(a)

odds = 0.37
print(odds / (1 + odds))

[1] 0.270073

(b)

p = 0.16
print(p / (1 - p))

[1] 0.1904762

Задача 10

Податоците од задачата се:

library(ISLR)

Weekly

(a)

Од нумеричката анализа забележуваме дека прирастот за 1, 2, 3, 4 и 5 недели веројатно има слична дистибуција (Lag1, Lag2, Lag3, Lag4, Lag5). Истото може да се забележи и за Today.

Поголем дел од прирастот е позитивен (просечниот и среден прираст за сите колони е над 0), што може да се објасни со тоа дека пазарот воглавно растел во 20те години кои се разгледани тука. Истото може и да се забележи со тоа што има повеќе Up од Down појавувања (55% се Up).

summary(Weekly)

      Year           Lag1               Lag2               Lag3         
 Min.   :1990   Min.   :-18.1950   Min.   :-18.1950   Min.   :-18.1950  
 1st Qu.:1995   1st Qu.: -1.1540   1st Qu.: -1.1540   1st Qu.: -1.1580  
 Median :2000   Median :  0.2410   Median :  0.2410   Median :  0.2410  
 Mean   :2000   Mean   :  0.1506   Mean   :  0.1511   Mean   :  0.1472  
 3rd Qu.:2005   3rd Qu.:  1.4050   3rd Qu.:  1.4090   3rd Qu.:  1.4090  
 Max.   :2010   Max.   : 12.0260   Max.   : 12.0260   Max.   : 12.0260  
      Lag4               Lag5              Volume            Today         
 Min.   :-18.1950   Min.   :-18.1950   Min.   :0.08747   Min.   :-18.1950  
 1st Qu.: -1.1580   1st Qu.: -1.1660   1st Qu.:0.33202   1st Qu.: -1.1540  
 Median :  0.2380   Median :  0.2340   Median :1.00268   Median :  0.2410  
 Mean   :  0.1458   Mean   :  0.1399   Mean   :1.57462   Mean   :  0.1499  
 3rd Qu.:  1.4090   3rd Qu.:  1.4050   3rd Qu.:2.05373   3rd Qu.:  1.4050  
 Max.   : 12.0260   Max.   : 12.0260   Max.   :9.32821   Max.   : 12.0260  
 Direction 
 Down:484  
 Up  :605  
           
           
           
           

Кога графички ќе ги претставиме податоците, забележуваме дека накај 950-от примерок се зголемува дисперзијата на прирастот. Ова е периодот околу 2008 година, кога се случи една од поголемите економски кризи.

plot(Weekly$Lag1, type="l", col="green", xlab="return", ylab="index")

plot(Weekly$Lag2, type="l", col="green", xlab="return", ylab="index")

plot(Weekly$Lag3, type="l", col="green", xlab="return", ylab="index")

plot(Weekly$Lag4, type="l", col="green", xlab="return", ylab="index")

plot(Weekly$Lag5, type="l", col="green", xlab="return", ylab="index")

(b)

Интерцептот е стратистички најзначаен. Неговата позитивна вредност укажува дека насоката има тенденција да е Up.

Lag2 е донекаде значајно, иако не премногу (сепак, p-вредноста не е претерано ниска). Ова укажува дека ако прирастот пред две недели бил позитивен, можеме да имаме некоја сигурност дека пазарот ќе се движи нагоре.

logistic_model = glm(
  Direction ∼ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
  data = Weekly,
  family = binomial 
)
summary(logistic_model)

Call:
glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
    Volume, family = binomial, data = Weekly)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6949  -1.2565   0.9913   1.0849   1.4579  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  0.26686    0.08593   3.106   0.0019 **
Lag1        -0.04127    0.02641  -1.563   0.1181   
Lag2         0.05844    0.02686   2.175   0.0296 * 
Lag3        -0.01606    0.02666  -0.602   0.5469   
Lag4        -0.02779    0.02646  -1.050   0.2937   
Lag5        -0.01447    0.02638  -0.549   0.5833   
Volume      -0.02274    0.03690  -0.616   0.5377   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1496.2  on 1088  degrees of freedom
Residual deviance: 1486.4  on 1082  degrees of freedom
AIC: 1500.4

Number of Fisher Scoring iterations: 4

(c)

Можеме да забележиме дека моделот инклинира повеќе кон предикции од типот Up. Од 484 примероци каде што насоката била Down, логистичкиот модел 430 пати направил false negative за Down класата и предвидел Up.

accuracy на моделот е 56.1%.

probas = predict(logistic_model, type = "response")
preds = rep("Down", length(probas))
preds[probas > 0.5] = "Up"
trues = Weekly$Direction
table(preds, trues)

      trues
preds  Down  Up
  Down   54  48
  Up    430 557

mean(preds == trues)

[1] 0.5610652

(d) (e) (f) (g) (h)

Ја правиме поделбата на тренирачко и тестирачко множество:

train_idxs = (Weekly$Year <= 2008)
train_weekly = Weekly[train_idxs, ]
test_weekly = Weekly[!train_idxs, ]
train_weekly

test_weekly

Можеме да ги направиме следниве заклучови за тестовите:

• логистичка регресија и LDA вадат исти резултати и се најдобри
• QDA учи скоро секогаш да вели Up, па е најлош од сите модели
• KNN со К=1 го има over-fittnato тренирачкото множество и не вади подобро од рандом

library(MASS)
library(class)
set.seed(1337)

get_conf_matrix <- function(name, preds, trues) {
  print(paste("Conf. Matrix for", name, sep=" "))
  print(table(preds, trues))
  print(paste("Accuracy for", name, sep=" "))
  mean(preds == trues)
}

get_logistic_conf_matrix <- function(name, model, data, trues) {
  probas = predict(model, data, type="response")
  preds = rep("Down", length(probas))
  preds[probas > 0.5] = "Up"
  get_conf_matrix(name, preds, trues)
}

get_da_conf_matrix <- function(name, model, data, trues) {
  preds = predict(model, data)$class
  get_conf_matrix(name, preds, trues)
}

get_knn_conf_matrix <- function(name, model, data, trues) {
  get_conf_matrix(name, model, trues)
}

trues = test_weekly$Direction

logistic_model = glm(Direction ∼ Lag2, data=train_weekly, family=binomial)
lda_model = lda(Direction ∼ Lag2, data=train_weekly)
qda_model = qda(Direction ∼ Lag2, data=train_weekly)
knn_model = knn(
  data.frame(train_weekly$Lag2), 
  data.frame(test_weekly$Lag2), 
  train_weekly$Direction, 
  k=1
)

get_logistic_conf_matrix(
  "logistic regression", 
  logistic_model, 
  test_weekly, 
  trues
)

[1] "Conf. Matrix for logistic regression"
      trues
preds  Down Up
  Down    9  5
  Up     34 56
[1] "Accuracy for logistic regression"
[1] 0.625

get_da_conf_matrix("LDA", lda_model, test_weekly, trues)

[1] "Conf. Matrix for LDA"
      trues
preds  Down Up
  Down    9  5
  Up     34 56
[1] "Accuracy for LDA"
[1] 0.625

get_da_conf_matrix("QDA", qda_model, test_weekly, trues)

[1] "Conf. Matrix for QDA"
      trues
preds  Down Up
  Down    0  0
  Up     43 61
[1] "Accuracy for QDA"
[1] 0.5865385

get_knn_conf_matrix("KNN 1", knn_model, test_weekly, trues)

[1] "Conf. Matrix for KNN 1"
      trues
preds  Down Up
  Down   21 29
  Up     22 32
[1] "Accuracy for KNN 1"
[1] 0.5096154

(i)

Ги испробуваме сите можни комбинации на предиктори (2^6). За секое подмножество испорбуваме логистичка регресија, LDA, QDA и 200 КНН класификатори. Вкупно се испробани 2^6 * (200 + 3) = 12992 модели во околу 2 минути. Од нив, највисока точност има KNN со К=122, а точноста изнесува 66.35%, користејќи го предикторот Lag1.

set.seed(1337)
predictors = c("Lag1", "Lag2", "Lag3", "Lag4", "Lag5", "Volume")

get_logistic_model_acc <- function(formula, train_data, test_data) {
  model = glm(as.formula(formula), data=train_data, family=binomial)
  probas = predict(model, test_data, type="response")
  preds = rep("Down", length(probas))
  preds[probas > 0.5] = "Up"
  return(mean(preds == trues))
}

get_da_model_acc <- function(formula, train_data, test_data, da) {
  model = da(as.formula(formula), data=train_data)
  preds = predict(model, test_data)$class
  return(mean(preds == trues))
}

get_knn_model_acc <- function(predictors, train_data, test_data, k) {
  preds = knn(
    train_data[predictors], 
    test_data[predictors], 
    train_data$Direction, 
    k=k
  )
  return(mean(preds == trues))
}

best_model_acc = 0
best_model_formula = ""
best_model_type = ""
best_model_k = ""

for (num_predictors in 1:length(predictors)) {
  predictor_combinations = combn(predictors, num_predictors)
  num_combinations = dim(predictor_combinations)[2]
  for (combination_idx in 1:num_combinations) {
    formula = paste(
      "Direction~",
      paste(predictor_combinations[, combination_idx], collapse="+"),
      sep=""
    )
    
    logistic_acc = get_logistic_model_acc(formula, train_weekly, test_weekly)
    if (logistic_acc > best_model_acc) {
      best_model_acc = logistic_acc
      best_model_formula = formula
      best_model_type = "logistic"
    }
    
    lda_acc = get_da_model_acc(formula, train_weekly, test_weekly, lda)
    if (lda_acc > best_model_acc) {
      best_model_acc = lda_acc
      best_model_formula = formula
      best_model_type = "LDA"
    }
    
    qda_acc = get_da_model_acc(formula, train_weekly, test_weekly, qda)
    if (qda_acc > best_model_acc) {
      best_model_acc = qda_acc
      best_model_formula = formula
      best_model_type = "QDA"
    }
    
    for (k in 1:200) {
      knn_acc = get_knn_model_acc(predictors, train_weekly, test_weekly, k)
      if (knn_acc > best_model_acc) {
        best_model_acc = knn_acc
        best_model_formula = formula
        best_model_type = "KNN"
        best_model_k = k
      }
    }
  }
}

print(paste("Best model: ", best_model_type))

[1] "Best model:  KNN"

print(paste("Acc: ", best_model_acc))

[1] "Acc:  0.663461538461538"

print(paste("Formula: ", best_model_formula))

[1] "Formula:  Direction~Lag1"

if (best_model_type == "KNN") {
  print(paste("K=", best_model_k))
}

[1] "K= 122"

Задача 11

Податоците од задачата се:

Auto

(a)

median_mpg = median(Auto$mpg)
mpg0 = rep(0, length(Auto$mpg))
mpg0[Auto$mpg > median_mpg] = 1
Auto$mpg0 = mpg0
Auto

(b)

Од boxplot-овите можеме да забележиме дека year и acceleration нема да бидат добри предиктори на mpg0. Останатите релативно добра сепарација прават помеѓу двете класи на mpg0. Најдобри се чинат weight, displacement и cylinders.

for (col in colnames(Auto)) {
  if (col == "mpg" || col == "mpg0" || col == "name") {
    next
  }
  
  boxplot(as.formula(paste(col, "~ mpg0")), data=Auto)
}

(c)

Правиме рандом поделба 80 - 20.

split = round(length(Auto$mpg) * 0.8)
idxs = sample(1:length(Auto$mpg))
train_idxs = idxs[1:split]
test_idxs = idxs[(split + 1):length(Auto$mpg)]
train_auto = Auto[train_idxs,]
test_auto = Auto[test_idxs,]
train_auto

test_auto

(d) (e) (f) (g)

set.seed(1337)

logistic_model = glm(
  mpg0 ~ weight + displacement + cylinders,
  data=train_auto, 
  family=binomial
)
probas = predict(logistic_model, test_auto, type="response")
preds = rep(0, length(probas))
preds[probas > 0.5] = 1
print(
  paste(
    "Logistic regression error rate:", 
    1 - mean(preds == test_auto$mpg0), 
    sep=" "
  )
)

[1] "Logistic regression error rate: 0.0897435897435898"

lda_model = lda(mpg0 ~ weight + displacement + cylinders, data=train_auto)
preds = predict(lda_model, test_auto)$class
print(paste("LDA error rate:", 1 - mean(preds == test_auto$mpg0), sep=" "))

[1] "LDA error rate: 0.0769230769230769"

qda_model = qda(mpg0 ~ weight + displacement + cylinders, data=train_auto)
preds = predict(qda_model, test_auto)$class
print(paste("QDA error rate:", 1 - mean(preds == test_auto$mpg0), sep=" "))

[1] "QDA error rate: 0.0641025641025641"

best_knn_k = 0
best_knn_acc = 0
for (k in 1:300) {
  preds = knn(
    train_auto[c("weight", "displacement", "cylinders")], 
    test_auto[c("weight", "displacement", "cylinders")], 
    train_auto$mpg0, 
    k=k
  )
  acc = mean(preds == test_auto$mpg0)
  if (acc > best_knn_acc) {
    best_knn_acc = acc
    best_knn_k = k
  }
}
print(paste(best_knn_k, "KNN error rate:", 1 - best_knn_acc, sep=" "))

[1] "5 KNN error rate: 0.0641025641025641"

Задача 12

Оваа задача е за едноставни функции и нема некаква поврзаност со материјата. Низ домашнава имам користено доста функции, така што ја скокнав оваа задача.

Задача 13

Ја додаваме медијаната:

median_crim = median(Boston$crim)
crim0 = rep(0, length(Boston$crim))
crim0[Boston$crim > median_crim] = 1
Boston$crim0 = crim0
Boston

Делиме на тренирачко и тестирачко множество:

split = round(length(Boston$crim0) * 0.8)
idxs = sample(1:length(Boston$crim0))
train_idxs = idxs[1:split]
test_idxs = idxs[(split + 1):length(Boston$crim0)]
train_boston = Boston[train_idxs,]
test_boston = Boston[test_idxs,]
train_boston

test_boston

Ги испробуваме сите можни комбинации на предиктори (2^13). За секое подмножество испробуваме логистичка регресија, LDA, QDA и 20 КНН класификатори. Вкупно се испробани 2^13 * (20 + 3) = 188416 модели во околу 7 минути. Од нив, највисока точност има логистичката регресија, а точноста изнесува 97.03%, користејќи ги предикторите indus, nox, rm, age, dis, rad, ptratio, black, lstat и medv.

Би испробале и повеќе опции за K, но времето на извршување расте многу (накај 3 часа би било за сите К од 1 до 200). За тоа нема потреба, бидејќи точноста и вака е доволно висока.

set.seed(1337)
predictors = c(
  "zn", "indus", "chas", "nox", "rm", "age", "dis", "rad", "tax", 
  "ptratio", "black", "lstat", "medv"
)
trues = test_boston$crim0

get_logistic_model_acc <- function(formula, train_data, test_data) {
  model = glm(as.formula(formula), data=train_data, family=binomial)
  probas = predict(model, test_data, type="response")
  preds = rep(0, length(probas))
  preds[probas > 0.5] = 1
  return(mean(preds == trues))
}

get_da_model_acc <- function(formula, train_data, test_data, da) {
  model = da(as.formula(formula), data=train_data)
  preds = predict(model, test_data)$class
  return(mean(preds == trues))
}

get_knn_model_acc <- function(predictors, train_data, test_data, k) {
  preds = knn(
    train_data[predictors], 
    test_data[predictors], 
    train_data$crim0, 
    k=k
  )
  return(mean(preds == trues))
}

best_model_acc = 0
best_model_formula = ""
best_model_type = ""
best_model_k = ""

for (num_predictors in 1:length(predictors)) {
  predictor_combinations = combn(predictors, num_predictors)
  num_combinations = dim(predictor_combinations)[2]
  for (combination_idx in 1:num_combinations) {
    formula = paste(
      "crim0~",
      paste(predictor_combinations[, combination_idx], collapse="+"),
      sep=""
    )
    
    logistic_acc = get_logistic_model_acc(formula, train_boston, test_boston)
    if (logistic_acc > best_model_acc) {
      best_model_acc = logistic_acc
      best_model_formula = formula
      best_model_type = "logistic"
    }
    
    lda_acc = get_da_model_acc(formula, train_boston, test_boston, lda)
    if (lda_acc > best_model_acc) {
      best_model_acc = lda_acc
      best_model_formula = formula
      best_model_type = "LDA"
    }

    qda_acc = get_da_model_acc(formula, train_boston, test_boston, qda)
    if (qda_acc > best_model_acc) {
      best_model_acc = qda_acc
      best_model_formula = formula
      best_model_type = "QDA"
    }

    for (k in 1:20) {
      knn_acc = get_knn_model_acc(predictors, train_boston, test_boston, k)
      if (knn_acc > best_model_acc) {
        best_model_acc = knn_acc
        best_model_formula = formula
        best_model_type = "KNN"
        best_model_k = k
      }
    }
  }
}

print(paste("Best model: ", best_model_type))

[1] "Best model:  logistic"

print(paste("Acc: ", best_model_acc))

[1] "Acc:  0.97029702970297"

print(paste("Formula: ", best_model_formula))

[1] "Formula:  crim0~indus+nox+rm+age+dis+rad+ptratio+black+lstat+medv"

if (best_model_type == "KNN") {
  print(paste("K:", best_model_k))
}