Resampling Methods

Задача 1

\[ Var(\alpha X+(1-\alpha)Y) =\\ \alpha^2Var(X)+(1-\alpha)^2Var(Y)+2\alpha(1-\alpha)Cov(X, Y)=\\ \alpha^2Var(X)+Var(Y)-2\alpha Var(Y)+\alpha^2Var(Y)+2\alpha Cov(X, Y) - 2\alpha^2Cov(X,Y) \\ \frac{\partial Var(\alpha X+(1-\alpha)Y)}{\partial \alpha} = \\ 2\alpha Var(X)-2Var(Y)+2\alpha Var(Y) + 2Cov(X, Y)-4\alpha Cov(X,Y)\\ \frac{\partial Var(\alpha X+(1-\alpha)Y)}{\partial \alpha} = 0 \implies\\ \alpha(Var(X)+Var(Y)-2Cov(X,Y))=Var(Y)-Cov(X,Y) \implies \\ \alpha=\frac{Var(Y)-Cov(X,Y)}{Var(X)+Var(Y)-2Cov(X,Y)}=\frac{\sigma_Y^2-\sigma_{XY}}{\sigma_X^2+\sigma_Y^2-2\sigma_{XY}} \]

Задача 2

(a)

Веројатноста дека првиот примерок во бутстрапот е ј-тиот примерок во оригиналот е \(\frac{1}{n}\). Веројатноста дека првиот примерок во бутстрапот не е ј-тиот примерок во оригиналот е \(1-\frac{1}{n}=\frac{n-1}{n}\).

(b)

Бидејќи влечеме со замена, второто влечење не зависи од првото, па веројатноста останува иста: \(1-\frac{1}{n}=\frac{n-1}{n}\).

(c)

Слично како и под б, за i-тиот примерок од бутстрапот важи дека веројатноста тој да не е ј-от примерок од оригиналот е \(1-\frac{1}{n}=\frac{n-1}{n}\). Сите влечења се независини, па затоа веројатноста да го нема ј-от примерок од оригиналот во целиот бутстрап е веројатноста сите примероци од бутстрапот да не се ј-от примерок од оригиналот што е \((1-\frac{1}{n})^n=(\frac{n-1}{n})^n\).

(d)

Веројатноста е обратна од тоа да го нема (под c), односно \(1-(1-\frac{1}{5})^5=1-\frac{4}{5}^5=0.67232\)

(e)

\(1-(1-\frac{1}{100})^{100}=1-\frac{99}{100}^{100}=0.6339677\)

(f)

\(1-(1-\frac{1}{10000})^{10000}=1-\frac{9999}{10000}^{10000}=0.632139\)

(g)

Конвергира кон 0.6321.

n=1:100000
p=1-((n-1)/n)^n
plot(n, p, type="l", ylim=c(0.5, 1))

(h)

Очекувано, резултатот е близок до 0.6321.

store = rep(NA, 10000)
for(i in 1:10000) {
  store[i] = sum(sample(1:100, rep=TRUE) == 4) > 0
}
mean(store)

[1] 0.6353

Задача 3

(a)

Тренирачкото множество се дели на к еднакви подмножества. Модел се тренира к пати. Во i-тото трениранје, i-тото подмножество се трга настрана и не се тренира на него, туку се пресметува точноста/error rate/. Error rate на тестирачкото множество се оценува како просек од error rates на секое од к-те крос-валидациски подмножества.

(b)

• cross validation има понизок bias од еден валидациски сет, но знае да биде поваријабилно
• LOOCV (cross validation со к=1) пак има најнизок bias, но затоа варијабилноста може да биде превисока, па може да сакаме да користиме cross validation; секако, LOOCV и доста поспоро се извршува

Задача 4

Користејќи го бутстрап методот, тренираме B модели. За секој модел, со бутстрап техниката, избираме n парови (X, Y) од постоечкото множество од n парови (X, Y) со замена. За секое вакво бутстрап подмножество, го тренираме моделот и правиме предикција за точното X за кое што сме заинтересирани. Добиваме оценувачи на Y \((Y_1^*, Y_2^*, ..., Y_B^*)\). Користејќи ја формулата (5.8) ја пресметуваме стандардната девијација.

Задача 5

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

library(ISLR)
Default

Ќе ја користиме оваа функција за да тренираме модел на тренирачко податочно множество и го валидираме на валидациско податочно множество, притоа враќајќи го error rate-от на валидациското податочно множество.

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

Ќе ја користиме оваа функција за да добиеме тренирачки и валидациски множества:

get_train_val_sets <- function(data, split_pct, num_val_sets) {
  set.seed(1337)
  idxs = 1:nrow(data)
  idxs_permutation = sample(idxs)
  val_size = round(nrow(data) * split_pct)
  train_sets = list()
  val_sets = list()
  for (val_idx in 1:num_val_sets) {
    val_seq = (1 + val_size * (val_idx - 1)):(1 + val_size * val_idx)
    train_sets[[val_idx]] = Default[idxs_permutation[-val_seq], ]
    val_sets[[val_idx]] = Default[idxs_permutation[val_seq], ]
  }
  return(list(train_sets, val_sets))
}

(a) (b)

sets = get_train_val_sets(Default, 0.2, 1)
train_set = sets[[1]][[1]]
val_set = sets[[2]][[1]]
get_logistic_model_err(
  "default ~ income + balance", 
  train_set, 
  val_set, 
  val_set$default
)

[1] 0.02398801

(c)

Резултатите од трите поделби се доста слични.

sets = get_train_val_sets(Default, 0.2, 3)
for (val_idx in 1:3) {
  train_set = sets[[1]][[val_idx]]
  val_set = sets[[2]][[val_idx]]
  print(
    get_logistic_model_err(
      "default ~ income + balance", 
      train_set, 
      val_set, 
      val_set$default
    )
  )
}

[1] 0.02398801
[1] 0.02398801
[1] 0.02448776

(d)

Нема некое подобрување.

sets = get_train_val_sets(Default, 0.2, 3)
for (val_idx in 1:3) {
  train_set = sets[[1]][[val_idx]]
  val_set = sets[[2]][[val_idx]]
  print(
    get_logistic_model_err(
      "default ~ income + balance + student", 
      train_set, 
      val_set, 
      val_set$default
    )
  )
}

[1] 0.02498751
[1] 0.02448776
[1] 0.02398801

Задача 6

(a)

model = glm(default ~ income + balance, family=binomial, data=Default)
coef(model)

  (Intercept)        income       balance 
-1.154047e+01  2.080898e-05  5.647103e-03 

(b)

set.seed(1337)
boot.fn <- function (data, idxs){
  return(
    coef(
      glm(default ~ income + balance, family=binomial, data=data[idxs, ])
    )
  )
}
boot.fn(Default, sample(nrow(Default), replace=TRUE))

  (Intercept)        income       balance 
-1.179284e+01  1.814249e-05  5.751864e-03 

(c)

library(boot)
set.seed(1337)
boot_stds = boot(Default, boot.fn, 100)
boot_stds

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Default, statistic = boot.fn, R = 100)


Bootstrap Statistics :
         original        bias     std. error
t1* -1.154047e+01 -4.699342e-02 4.513118e-01
t2*  2.080898e-05  3.373744e-07 5.121736e-06
t3*  5.647103e-03  1.530472e-05 2.237195e-04

(d)

Ги плотираме коефициентите добиени со glm (зелена боја) и отстапувањата (портокалова боја) заедно со просекот од бутстрапот (сина боја, испрекината линија) и отстапувањата од бутстрапот (црвена боја, испрекината). Се забележува дека вредностите се доста слични.

for (coef_idx in 1:3) {
  min_est = min(boot_stds$t[, coef_idx])
  max_est = max(boot_stds$t[, coef_idx])
  range = max_est - min_est
  padding = range * 0.1
  ylim = c(min_est - padding, max_est + padding + range / 3)
  plot(
    boot_stds$t[, coef_idx], 
    type="l", 
    col="black", 
    ylab=paste("beta", coef_idx, sep=""),
    xlab="Bootstrap iteration",
    ylim=ylim
  )
  abline(a=mean(boot_stds$t[, coef_idx]), b=0, col="blue", lty=2)
  abline(
    a=mean(boot_stds$t[, coef_idx]) + sd(boot_stds$t[, coef_idx]), 
    b=0, 
    col="red", 
    lty=2
  )
  abline(
    a=mean(boot_stds$t[, coef_idx]) - sd(boot_stds$t[, coef_idx]), 
    b=0, 
    col="red", 
    lty=2
  )
  abline(a=coef(model)[coef_idx], b=0, col="green")
  model_std = summary(model)$coefficients[coef_idx, "Std. Error"]
  abline(a=coef(model)[coef_idx] + model_std, b=0, col="orange")
  abline(a=coef(model)[coef_idx] - model_std, b=0, col="orange")
  legend(
    0, 
    ylim[2], 
    legend=c("Model Beta", "Model Std"), 
    col=c("green", "orange"), 
    lty=c(1, 1)
  )
  legend(
    40, 
    ylim[2], 
    legend=c("Mean Beta", "Std Beta"), 
    col=c("blue", "red"), 
    lty=c(2, 2)
  )
}

Задача 7

(a)

full_model = glm(Direction ~ Lag1 + Lag2, family=binomial, data=Weekly)
summary(full_model)

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

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.623  -1.261   1.001   1.083   1.506  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.22122    0.06147   3.599 0.000319 ***
Lag1        -0.03872    0.02622  -1.477 0.139672    
Lag2         0.06025    0.02655   2.270 0.023232 *  
---
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: 1488.2  on 1086  degrees of freedom
AIC: 1494.2

Number of Fisher Scoring iterations: 4

(b)

all_but_first = glm(
  Direction ~ Lag1 + Lag2, 
  family=binomial, 
  data=Weekly[2:nrow(Weekly), ]
)
summary(all_but_first)

Call:
glm(formula = Direction ~ Lag1 + Lag2, family = binomial, data = Weekly[2:nrow(Weekly), 
    ])

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6258  -1.2617   0.9999   1.0819   1.5071  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.22324    0.06150   3.630 0.000283 ***
Lag1        -0.03843    0.02622  -1.466 0.142683    
Lag2         0.06085    0.02656   2.291 0.021971 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1494.6  on 1087  degrees of freedom
Residual deviance: 1486.5  on 1085  degrees of freedom
AIC: 1492.5

Number of Fisher Scoring iterations: 4

(c)

Не е точно класифицирана.

proba_first = predict(all_but_first, Weekly[1:1, ], type="response")
pred = "Down"
if (proba_first > .5) {
  pred = "Up"
}
proba_first

        1 
0.5713923 

Weekly[1, "Direction"] == pred

[1] FALSE

(d) (e)

Оценкето е дека 44% од предвидувањата се погрешни. Сепак, овој пристап дава оценка со висока варијабилност, па не можеме да сме комплетно сигурни дека моделот е толку лош.

loocv_err = 0

for (idx in 1:nrow(Weekly)) {
  model = glm(Direction ~ Lag1 + Lag2, family="binomial", data=Weekly[-idx, ])
  proba = predict(model, Weekly[idx, ], type="response")
  pred = "Down"
  if (proba > 0.5) {
    pred = "Up"
  }
  if (pred != Weekly[idx, "Direction"]) {
    loocv_err = loocv_err + 1 / nrow(Weekly)
  }
}

print(loocv_err)

[1] 0.4499541

Задача 8

(a)

n (бројот на примероци) e 100, a p (бројот на предиктори) е 1. Формулата е: \(y =-2x^2+x+err\).

set.seed(1)
y = rnorm(100)
x = rnorm(100)
y = x - 2 * x ^ 2 + rnorm(100)

(b)

Делува како да има квадратна врска помеѓу x и y (реално, знаеме дека има, од моделот).

plot(x, y)

(c) (d)

Резултатите се очекувано исти. Бидејќи користиме LOOCV, секој можен датасет од 99 елементи го пробуваме еднаш. Нема простор за случајност тука, па резултатите не зависат од seed-от.

data = data.frame(x, x ^ 2, x ^ 3, x ^ 4, y)

loocv_sint <- function(seed) {
  set.seed(seed)
  
  model = glm(y ∼ x, data=data)
  loocv_err = cv.glm(data, model)
  print(loocv_err$delta)
  
  model = glm(y ∼ x + x.2, data=data)
  loocv_err = cv.glm(data, model)
  print(loocv_err$delta)
  
  model = glm(y ∼ x + x.2 + x.3, data=data)
  loocv_err = cv.glm(data, model)
  print(loocv_err$delta)

  model = glm(y ∼ x + x.2 + x.3 + x.4, data=data)
  loocv_err = cv.glm(data, model)
  print(loocv_err$delta)
}

for (seed in 1:2) {
  loocv_sint(seed)
  print("----------")
}

[1] 5.890979 5.888812
[1] 1.086596 1.086326
[1] 1.102585 1.102227
[1] 1.114772 1.114334
[1] "----------"
[1] 5.890979 5.888812
[1] 1.086596 1.086326
[1] 1.102585 1.102227
[1] 1.114772 1.114334
[1] "----------"

(e)

Квадратниот има најниско LOOCV. Тоа е очекувано, бидејќи оригиналниот модел со кој е добиено податочното множество е квадратен.

(f)

Во секој модел, освен линеарниот, најважни се \(x\) и \(x^2\). Тоа е очекувано, бидејќи податоците се добиени со квадратен модел. Линеарниот модел, нормално, доста лошо ги учи овие квадратни податоци, па кај него само интерцептот е важен.

summary(glm(y ∼ x, data=data))
summary(glm(y ∼ x + x.2, data=data))
summary(glm(y ∼ x + x.2 + x.3, data=data))
summary(glm(y ∼ x + x.2 + x.3 + x.4, data=data))

Задача 9

Оваа задача има се од претходните во неа, па затоа е скокната.