Motivation

Given a prediction function \(\hat{f}\), that we trained on a training dataset \(\mathcal{D}_n = \{(X_1, Y_1), \ldots, (X_n, Y_n)\}\), we would like to:

• Evaluate its performance, i.e. its risk \[ R(\hat{f}) = \mathrm{E}_{\mathcal{D}_n, Y, X}[\ell(Y, \hat{f}(X))] \] loss functions typically \(\ell(y,f) = (y-f)^2\) (regression task), \(\ell(y,f) = \mathbb{1}\{y \neq f\}\) (classification task)

• Choose the hyperparameters that minimise test error:

  • \(\hat{f}(x)\) often depends on hyperparameters, i.e. parameters that we do not choose in the training phase of \(\hat{f}\).
  • Examples: the number of neighbours in kNN, the learning rate \(\alpha\) in gradient descent for logistic classification

Problem

We only have one dataset \(\mathcal{D}_n\), which we already used for training \(\hat{f}\).

We can look at the training error:

\[ \mathrm{Training \, error} = \frac{1}{n} \sum_{i = 1}^n \ell(Y_i, \hat{f}(X_i)) \,. \]

But, because we used the training data already for fitting \(\hat{f}\), our estimate will be too optimistic. If \(\hat{f}\) is powerful enough, it can achieve zero training error.

Validation set approach

We can split our data \(\mathcal{D}_n\) into two sets, a training set and a validation set.

We train \(\hat{f}\) on only set of the data, and evaluate on the other.

1. Randomly split data into train and test set
2. Train \(\hat{f}\) on training set
3. Estimate test error on test set

Validation set approach

• Random split introduces uncertainty

• Less training data → worse training → worse predictions → overestimate test error

Validation set approach: Model selection

Say our prediction \(\hat{f}\) has a hyperparameter, call it \(\lambda\), that we want to tune.

We do the following:

1. Split data into train and validation set
2. For a range of hyperparmeters \(\lambda_i \in \{\lambda_1, \ldots, \lambda_k\}\), do:
   1. Train \(\hat{f}(\cdot;\lambda_i)\) on training data
   2. Evluate error on validation data
3. Pick the model \(\hat{f}(\cdot; \hat{\lambda})\), where \(\hat{\lambda}\) achieves smallest error on validation set

Is the error on the validation set a good estimate for the test error of \(\hat{f}(\cdot; \hat{\lambda})\)?

Illustration: Auto data

library(ISLR) 

str(Auto)

'data.frame':   392 obs. of  9 variables:
 $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
 $ cylinders   : num  8 8 8 8 8 8 8 8 8 8 ...
 $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
 $ horsepower  : num  130 165 150 150 140 198 220 215 225 190 ...
 $ weight      : num  3504 3693 3436 3433 3449 ...
 $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
 $ year        : num  70 70 70 70 70 70 70 70 70 70 ...
 $ origin      : num  1 1 1 1 1 1 1 1 1 1 ...
 $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...

plot(mpg ~ horsepower, data = Auto)

Illustration: Auto data

We want to fit a model predicting mpg as a function of horsepower

Prediction function is a polynomial of degree \(d\) of horsepower. I.e.

• \(d = 1\): \(\hat{f}(x; 1) = \beta_0 + \beta_1 x\)
  
• \(d = 2\): \(\hat{f}(x; 2) = \beta_0 + \beta_1 x + \beta_2 x^2\)
• \(d = 3\): \(\hat{f}(x; 2) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3\)
• …

Use validation set approach to pick best hyperparameter and report estimate of test error under squared loss

Illustration: Auto data

Train, validate, test split

• Train: half
• Validate, test: one quarter

set.seed(2022)
n <- nrow(Auto)
shuffle <- sample(n) 
Auto_train <- Auto[shuffle[1:(n/2)], ]
Auto_validate <- Auto[shuffle[(n/2+1):(3*n/4)], ]
Auto_test <- Auto[shuffle[(3*n/4+1):n], ]
dim(Auto_train)

[1] 196   9

Illustration: Auto data

Hyperparameter grid

\[ d_i \in \{1, 2, \ldots, 10\}\]

MSE <- function(u, v){
mean((u - v)^2)
}

max_deg <- 10
MSE_train <- rep(0, max_deg) 
MSE_validate <- rep(0, max_deg) 
for (deg in 1:max_deg){
Auto_poly = lm(mpg ~ poly(horsepower, deg), data=Auto_train)

pred_train <- predict(Auto_poly, newdata=Auto_train)
pred_validate <- predict(Auto_poly, newdata=Auto_validate)

MSE_train[deg] <- MSE(pred_train, Auto_train$mpg)
MSE_validate[deg] <- MSE(pred_validate, Auto_validate$mpg)
}
plot(c(0,max_deg), range(c(MSE_train, MSE_validate)), type='n', 
   xlab='deg', ylab='MSE')
points(1:max_deg, MSE_train, col='blue', type='b')
points(1:max_deg, MSE_validate, col='orange', type='b')
legend("topright", lty = c(1,1), col = c("blue", "orange"), legend = c("Train", "Validate"))

Illustration: Auto data

Report estimated test error

best_deg <- which.min(MSE_validate)

Auto_poly <-  lm(mpg ~ poly(horsepower, best_deg), data=Auto_train)
pred_test <- predict(Auto_poly, newdata=Auto_test)
best_mse <- MSE(pred_test, Auto_test$mpg)

cat("Test error estimate for d = ", best_deg, ": ", best_mse)

Test error estimate for d =  4 :  18.93462

Activity

Repeat the same validation set approach on another dataset, e.g. the house data from the MASS package

library(MASS) 
str(Boston)

'data.frame':   506 obs. of  14 variables:
 $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
 $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
 $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
 $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
 $ rm     : num  6.58 6.42 7.18 7 7.15 ...
 $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
 $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
 $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
 $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
 $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
 $ black  : num  397 397 393 395 397 ...
 $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
 $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Leave-one-out cross validation

Aim to mitigate the variability of the validation set approach and be more efficient with our data

\[ \text{LOOCV}(\hat{f}) = \frac{1}{n} \sum_{i = 1}^n \ell(Y_i, \hat{f}^{(-i)}(X_i)) \]

• Each split like train, validation split

Leave-one-out cross validation

• Low bias: Beacuse train on almost foll data, LOOCV is a good (low bias) approximation of the error of \(\hat{f}\), i.e.  \[ \mathrm{E}_{\mathcal{D}_n}[\text{LOOCV}(\hat{f})] \approx \mathrm{E}_{X,Y,\mathcal{D}_n}[\ell(Y, \hat{f}(X))] \]

• High variance:

  • Because each summands predicts only a single point, we get a noisy estimate of expected loss
  • Because training sets are very similar, high correlation between them

  → If datasets change, we expect to see large fluctuations in LOOCV

• Computationally intensive: We need to fit a model for \(n\) different data folds

K-fold cross validation

Compromise between validation set approach and LOOCV:

1. Partition data \(\mathcal{D}_n\) into \(K\) cells \(C_1, \ldots, C_K\)
2. For fold \(k\), train \(\hat{f}^{(-k)}\) on \(\mathcal{D}_n \setminus C_k\), validate on \(C_k\)
3. Aggregate: \[ \text{CV}_k(\hat{f}) = \frac{1}{n} \sum_{k = 1}^K \sum_{j \in C_i} \ell(Y_j, \hat{f}^{(-k)}(X_j)) \]

K-fold cross validation

• Split uncertainty: K-fold CV still uses random data split, but variability is substantially lower than for validation set: Overlap between folds reduces impact of split

• Bias: Typically less bias than validation set approach, more bias than LOOCV
• Variance: Less variance than both LOOCV, and validation set approach

→ Empirically: K-fold CV with \(5-10\) folds: Sweet spot between variance, bias, and computational effort

Illustration: Auto data

Use library boot for CV

library(boot) #for cv.glm
Auto_train <- rbind(Auto_train, Auto_validate)
dim(Auto_train)

[1] 294   9

Illustration: Auto data

Run LOOCV, 5-fold CV, and 10-fold CV to estimate

\[ \hat{f}(x;d) = \beta_0 + \sum_{i = 1}^d \beta_i x^i \]

num <- 5
cv_error1 = rep(0, num)
cv_error5 = rep(0, num)
cv_error10 = rep(0, num)
degree = 1:num
for (deg in degree) {
glm_fit <-  glm(mpg ~ poly(horsepower, deg), data = Auto)
## CV errors using naive LOOCV
cv_error1[deg] = cv.glm(Auto, glm_fit)$delta[1]
## CV errors using K-fold-CV
cv_error5[deg] = cv.glm(Auto, glm_fit, K = 5)$delta[1]
## CV errors using K-fold-CV
cv_error10[deg] = cv.glm(Auto, glm_fit, K = 10)$delta[1]
}

plot(degree, cv_error1, type = "b", col = "black")
lines(degree, cv_error5, type = "b", col = "red")
lines(degree, cv_error10, type = "b", col = "blue")
legend("topright", legend = c("LOOCV", "5-fold CV", "10-fold CV"), col = c("black", "red", "blue"), lwd = 2)

Activity

Repeat the above with a different set of predictors.

Bootstrap

Suppose we are given data \((Y_1,X_1), \ldots, (Y_n, X_n)\) that is an i.i.d. sample from some unknown distribution \(P\).

We compute an estimate \(\hat \theta\) of some quantity of interest, e.g. the MSE of a ML-algorithm. We would like to know something about its distribution

Monte carlo method

If we had access to \(P\), we could:

1. Draw \(B\) independent training samples \(\mathcal{D}_1,\ldots, \mathcal{D}_B\)
2. Get estimates \(\hat{\theta}_1,\ldots, \hat{\theta}_B\)
3. Use the empirical distribution of our estimates as an approximation of the distribution of the estimator itself

Bootstrap

Empirical distribution

The empirical distrubtion function of a sample of \(n\) i.i.d. random variables \(X\) is

\[ P_n(x) = \frac{1}{n} \sum_{i = 1}^n \mathbb{1}\{X_i \leq x\} \]

set.seed(123) 
N <- 50

xs <- sort(rnorm(N)) 

emp_cdf <- 1 : N / N

x_grid <- seq(from = -3, to = 3, length.out = 1000)

plot(xs, emp_cdf, col = "blue", lwd = 2, type = "l", main = "Empirical CDF vs. CDF") 
lines(x_grid, pnorm(x_grid), col = "green", lwd = 2) 
legend("bottomright", col = c("blue", "green"), lwd = rep(2,2), legend = c("Emp. CDF", "CDF")) 

Bootstrap

Bootstrap idea

Rather than using the unknown distribution \(P\) to sample data, we use the empirical distribution \(P_n\) to sample data.

• This corresponds to sampling from our training data \(\mathcal{D}\) uniformly at random with replacement

Illustration: Portfolio data

The Portfolio data set is in the ISLR package described in Section 5.2.

'data.frame':   100 obs. of  2 variables:
 $ X: num  -0.895 -1.562 -0.417 1.044 -0.316 ...
 $ Y: num  -0.235 -0.885 0.272 -0.734 0.842 ...

Illustration: Portfolio data

Wish to invest a fixed sum of money in two financial assets that yield returns of \(X\) and \(Y\), respectively, where X and Y are random quantities. We will invest a fraction \(\alpha\) of our money in \(X\), and will invest the remaining \(1 - \alpha\) in \(Y\).

The formula for \(\alpha\) that results in the minimum risk investment of form \(\alpha X + (1-\alpha)Y\) is

\[ \alpha = \frac{\text{Var}(Y) - \text{Cov}(X,Y)}{\text{Var}(X)+ \text{Var}(Y)- 2\text{Cov}(X,Y)} \]

Illustration: Portfolio data

To illustrate the use of the bootstrap on this data, we must first create:

• a function, alpha_fn() that takes as input the \((X,Y)\) data and returns an estimate of alpha using the estimator (5.7) to the observations corresponding to the rows of the data index.
• a vector indicating which observations should be used to estimate \(\alpha\). Note that \(\alpha = \arg \min \text{var}[\alpha X + (1 - \alpha Y)]\).

alpha_fn <-  function(data, index) {
X <-  data$X[index]
Y <-  data$Y[index]
# solution of alpha s.t. min. var
alpha <-  (var(Y) - cov(X, Y)) / (var(X) + var(Y) - 2 * cov(X, Y))
return(alpha)
}

randomly select \(100\) observations from \(1\) to \(100\) with replacement, i.e. construct a new bootstrap data and compute the corresponding alpha

set.seed(1)
alpha_fn(data <-  Portfolio, index = sample(1:100, size = 100, replace = T))

[1] 0.7368375

Illustration: Portfolio data

We can implement a bootstrap analysis by performing this command many times, recording all of the corresponding estimates for α, and computing the resulting standard. deviation.

Bootstrap using for loop:

B <- 1000
boot_alpha = rep(0, B)
for(i in 1:B){
boot_index <-  sample(x = 1:100, size = 100, replace = T)
boot_alpha[i] <-  alpha_fn(data = Portfolio, index = boot_index)
}
mean(boot_alpha)

[1] 0.5739061

sd(boot_alpha)

[1] 0.09285792

boot automates this approach. To produce 1000 bootstrap estimates for alpha

boot(Portfolio, alpha_fn, R = 1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Portfolio, statistic = alpha_fn, R = 1000)


Bootstrap Statistics :
     original       bias    std. error
t1* 0.5758321 -0.001539618  0.08846251

Illustration: Auto data

Now we will find the bootstrap estimates of the standar error of the linear regression coefficients from the model of mpg on horsepower in Auto data

We get the estimated standard errors from lm:

Call:
lm(formula = mpg ~ horsepower, data = Auto)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.5710  -3.2592  -0.3435   2.7630  16.9240 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 39.935861   0.717499   55.66   <2e-16 ***
horsepower  -0.157845   0.006446  -24.49   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.906 on 390 degrees of freedom
Multiple R-squared:  0.6059,    Adjusted R-squared:  0.6049 
F-statistic: 599.7 on 1 and 390 DF,  p-value: < 2.2e-16

Now we will use bootstrap to estimate the standard errors of the coefficients:

Illustration: Auto data

First define a function that returns the intercept and slope estimates for the linear regression model:

(Intercept)  horsepower 
 39.9358610  -0.1578447 

(Intercept)  horsepower 
 40.9483172  -0.1636591 

Now we use boot() to compute the standard errors of 1000 bootstrap estimates for the intercept and slope

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot_fn, R = 1000)


Bootstrap Statistics :
      original        bias    std. error
t1* 39.9358610  0.0141715259 0.842049552
t2* -0.1578447 -0.0002181104 0.007203328

Illustration: Auto data

Compare with standard formula results for the regression coefficients in a linear model

              Estimate  Std. Error   t value      Pr(>|t|)
(Intercept) 39.9358610 0.717498656  55.65984 1.220362e-187
horsepower  -0.1578447 0.006445501 -24.48914  7.031989e-81

What can you conclude from the different results?

There is some difference between the SE estimates of the parameters using bootstrap approach and the formulas coming from assuming a linear model. Standard error estimates depend on \(\sigma^2\) that first needs to be estimated. \(\hat{\sigma}^2\) depends on the assumption that a linear model is correct, if it is not, then the residuals from the linear fit will be inflated, and consequently so will \(\hat{\sigma}^2\). The bootstrap approach does not rely on these assumptions so might provide more accurate estimates.

Activity

Redo everything for polynomial regression with degree=2

Check your understanding

In cross-validation, what is the main purpose of repeatedly splitting the data into training and validation sets?

Summary

Cross-validation

• Used to estimate test-errors
• \(K\)-fold CV for \(K=5\) to \(K=10\) sweetspot between bias, variance, and computational effort

Bootstrap

• Used to estimate distributional properties of an estimator whose distribution is unknown
• Use to construct estimates of standard errors, CIs

Key concepts

• Linear classifiers, LDA, logistic classifier (Q1 a, Q2)
• Loss, surrogate loss (Q1 b, c)
• Discriminant based classifiers (Q2, Q3, Q4 a)
  
• ROC (Q4 b)

Question 1

In this exercise, we consider binary classification, where \(Y = −1\) denotes a negative example and \(Y = 1\) denotes a positive example.

1. Show that the output of a logistic regression classifier with parameter \(\hat{\beta}\) can be expressed as follows \[ \psi(x) = H(x^\top \hat\beta) \,, \] where \(H(z) = 1\) if \(z \geq 0\) and \(H(z) = 0\) if \(z < 0\).

Question 1a

H <- function(z){
    ifelse(z < 0, 0, 1) 
}

zs <- seq(from = -5, to = 5, length.out= 1000) 

plot(zs, H(zs), col = "blue", lwd = 2, type = "l", main = expression(H(z)))

plot(zs, plogis(zs), col = "red", lwd = 2, type = "l", main = "Logistic classifier")
abline(h = 0.5, col = "grey", lty = 2, lwd = 2) 
lines(zs, H(zs), col = "blue", lwd = 2)

Question 1a

Show

Logistic classifier is a \(0\)-\(1\) step function in \(x^\top\hat{\beta}\)

Plug-in

\[\begin{equation} \begin{aligned} \psi(x) &= \mathbb{1}\{ \hat{\eta}(x) \geq 1 / 2 \} \\ &= \mathbb{1}\{ 1 / (1 + \exp\{-x^\top \hat \beta \}) \geq 1 / 2 \} \\ &= \mathbb{1}\{ 2 \geq 1 + \exp\{-x^\top \hat \beta \} \} \\ &= \mathbb{1}\{ 1 \geq \exp\{-x^\top \hat \beta \} \} \\ &= \mathbb{1}\{ 0 \geq -x^\top \hat \beta \} \\ &= \mathbb{1}\{ x^\top \hat \beta \geq 0 \} \\ &= H(x^\top \hat \beta) \,. \end{aligned} \end{equation}\]

Learnings

• Go from definition to statement

Question 1b

Let \(\ell: \Re \to \Re\) with \(\ell(z) = \log(1 + \exp\{-z\})\). Show that for logistic regression, the negative log-likelihood (referred to as binary cross-entropy or log loss) is \(\ell(-yx^\top \beta)\).

Question 1b

log-likelihhod:

\[ \mathcal{L} = \mathbb{1}\{Y = 1\} \log(\Pr(Y = 1 \mid X = x)) + \mathbb{1}\{Y = -1\}\log(\Pr(Y = -1 \mid X = x)) \]

\(y = 1\)

\[ \begin{aligned} \mathcal{L} & = \log(\Pr(Y = 1 \mid X = x))\\ &= \log\left(\frac{1}{1 + \exp\{-x^\top \hat \beta\}} \right) \\ &= \log(1) - \log \left({1 + \exp\{-x^\top \hat \beta\}} \right) \\ &= - \ell(-y x^\top \beta) \end{aligned} \]

Question 1b

\(y = -1\)

\[ \begin{aligned} \mathcal{L} & = \log(\Pr(Y = 1 \mid X = x))\\ &= \log\left(\frac{\exp\{-x^\top \beta\}}{1 + \exp\{-x^\top \hat \beta\}} \right) \\ &= \log\left(\frac{1}{1 + 1 / \exp\{-x^\top \hat \beta\}} \right) \\ &= \log\left(\frac{1}{1 + \exp\{x^\top \hat \beta\}} \right) \\ &= - \log(1 + \exp\{x^\top \hat \beta\}) \\ &= - \ell(-y x^\top \beta) \end{aligned} \]

Question 1b

Learnings

• Go from definition to statement
• log-likelihood estimation minimises cross-entropy loss

Question 1c

Consider a linear classifier \(H(x^\top \hat \beta)\) with the loss \(\ell(yx^\top \hat{\beta})\), using either the \(0\)-\(1\) loss function \(\ell(z) = \mathbb{1}\{z \leq 0\}\), the hinge loss function \(\ell(z) = \max \{0, 1-z\}\), or the squared hinge loss function \(\ell(z) = \max \{0, 1-z\}^2\).

1. Compare the values of the three loss functions over the domain R and discuss their continuity properties.
2. Which of the three loss functions is differentiable at each point \(z \in \Re\)?
3. How does the binary cross-entropy function compare to the above three functions (up to scaling with a multiplicative constant)?

Question 1c

hinge <- function(z){
    pmax(0, 1 - z) 
}

sq_hinge <- function(z){
    pmax(0, 1 - z)^2
}

z_o <- function(z){
    ifelse(z <= 0, 1, 0) 
}

plot(zs, hinge(zs), col = "black", lwd = 2, type = "l") 
lines(zs, sq_hinge(zs), col = "blue", lwd = 2)  
lines(zs, z_o(zs), col = "red", lwd = 2)  
legend("topright", col = c("black", "blue", "red"), lwd = rep(2,3), legend = c("hinge", "sq. hinge", "0-1")) 

• Hinge and squared hinge upper bound the \(0\)-\(1\) loss
• \(0\)-\(1\) loss has discontinuity at \(0\), others continuous
• Squared hinge loss is differentiable on \(\Re\), \(0\)-\(1\) loss is not differentiable at \(0\), hinge loss not differentiable at \(1\)

Question 1c

cross_entropy <- function(z, a = 1 / log(2)){
    a * log(1 + exp(-z)) 
}


plot(zs, hinge(zs), col = "black", lwd = 2, type = "l") 
lines(zs, sq_hinge(zs), col = "blue", lwd = 2)  
lines(zs, z_o(zs), col = "red", lwd = 2)  
lines(zs, cross_entropy(zs), col = "orange", lwd = 2)
legend("topright", col = c("black", "blue", "red", "orange"), lwd = rep(2,43), legend = c("hinge", "sq. hinge", "0-1", "cross entropy")) 

• upper bounds \(0\)-\(1\) loss for \(a > 1 / \log(2)\)
• can be seen as smooth approximation to hinge-loss function

Question 1c

Learnings:

• understand properties of difference loss functions
• maximum-likelihood estimate in logistic regression can be seen as smooth approximation to minimising empirical \(0\)-\(1\) loss (surrogate loss)

Question 2a

Consider LDA for \(p = 1\). a. Show that the discriminator can be expressed as follows: \[ \delta_k(x) = x \frac{\hat{\mu}_k}{\hat{\sigma}^2} - \frac{\hat{\mu}_k^2}{2 \hat{\sigma}^2} + \log(\hat{\pi}_k) \]

Recap

Discriminant based classifiers

Classifier \(\psi: \Re^p \to \{1 , \ldots, K\}\) is a discriminant based classifier if it can be expressed as

\[ \psi(x) = \underset{k}{\arg \max} \, \delta_k(x) \,, \] where \(\delta_k: \Re^p \to \Re\) is called the discriminant function. Any monotone transformation of \(\delta_k\) gives the same predictions!

Bayes discriminant function under \(0\)-\(1\) loss

\[ \eta_k(x) = \Pr(Y = k \mid X = x) = \frac{\pi_k f_k(x)}{\sum_{k = 1}^K \pi_k f_k(x)} \]

• Bayes’ Theorem: \(\Pr(A \mid B) = \frac{\Pr(A) \Pr(B \mid A)}{\Pr(B)}\)
• Law of total probablity: \(\Pr(B) = \sum_{i = 1}^n \Pr(A) \Pr(B_n \mid A)\), for \(B_n\) disjoint with \(\bigcup_{i =1}^n B_n = B\)

Recap

LDA

Assumes that \(X \mid Y = k\) is \(\mathrm{N}(\mu_k, \Sigma)\), i.e. 

\[ f_k(x) = \frac{1}{(2\pi)^{p/2}(\det\Sigma)^{1/2}} e^{-\frac{1}{2}(x-\mu_k)^\top\Sigma^{-1}(x-\mu_k)} \]

library(MASS)
library(ggplot2)
library(ggnewscale)

set.seed(42)

mu1 <- c(2, 3)
mu2 <- c(7, 8)
Sigma <- matrix(c(2, 1, 1, 2), 2, 2)
n <- 1000

X1 <- as.data.frame(mvrnorm(n, mu1, Sigma))
X2 <- as.data.frame(mvrnorm(n, mu2, Sigma))
colnames(X1) <- colnames(X2) <- c("x", "y")

X1s <- X1[sample(1:n, n/10), ]
X2s <- X2[sample(1:n, n/10), ]

ggplot() +
  stat_density_2d(
    data = X1, aes(x = x, y = y, fill = after_stat(level)),
    geom = "polygon", color = NA, alpha = 0.5
  ) +
  scale_fill_gradient(low = "#c6dbef", high = "#08519c", guide = "none") +
  
  new_scale_fill() +
  
  stat_density_2d(
    data = X2, aes(x = x, y = y, fill = after_stat(level)),
    geom = "polygon", color = NA, alpha = 0.5
  ) +
  scale_fill_gradient(low = "#c7e9c0", high = "#006d2c", guide = "none") +

  geom_point(data = X1s, aes(x = x, y = y), color = "#08519c", alpha = 0.8, size = 0.9) +
  geom_point(data = X2s, aes(x = x, y = y), color = "#006d2c", alpha = 0.8, size = 0.9) +
  
  theme_minimal(base_size = 14) +
  theme(legend.position = "none") +
  labs(
    title = expression("LDA modelling assumption"),
    x = expression(X[1]),
    y = expression(X[2])
  )

Recap LDA

With Bayes discriminant function

\[ \log (\pi_k f_k(x)) = \log\pi_k - \frac{1}{2}(x-\mu_k)^\top \Sigma^{-1}(x-\mu_k) \]

LDA estimates \(\pi_k, \mu_k, \Sigma\) to get

\[ \delta_k^{\mathrm{LDA}}(x):= \log\hat \pi_k - \frac{1}{2}(x-\hat \mu_k)^\top \hat\Sigma^{-1}(x-\hat \mu_k) \]

The learnable parameters \(\pi_k, \Sigma, \mu_k\) are estimated from the data:

• \(\hat{\pi}_k = \frac{n_k}{n} = \frac{1}{n} \sum_{i = 1}^n \mathbb{1}\{Y_i = k\}\)
• \(\hat\mu_k = \frac{1}{n_k}\sum_{i = 1}^n \mathbb{1}\{Y_i = k\} X_i\)
• \(\hat{\Sigma} = \frac{1}{n - K} \sum_{k = 1}^K \sum_{i = 1}^n \mathbb{1}\{Y_i = k\} (X_i - \hat{\mu}_k)(X_i - \hat{\mu}_k)^\top\)

Question 2a

Show

For \(p = 1\)

\[ \delta_k(x) = x \frac{\hat{\mu}_k}{\hat{\sigma}^2} - \frac{\hat{\mu}_k^2}{2 \hat{\sigma}^2} + \log(\hat{\pi}_k) \]

Know

• \(\delta_k(x) = \log(\hat\pi_k \hat f_k(x)) + c\), where \(c\) does not depend on class
  
• \(\hat{f}_k(x) = \frac{1}{\sqrt{2 \pi \hat{\sigma}^2}} \exp\{-(x - \mu_k)^2/ 2 \hat\sigma^2 \}\)

Question 2a

Plug in

\[ \begin{aligned} \log(\hat\pi_k \hat f_k(x)) &=\log(\hat{\pi}_k) -\frac{1}{2}\log(2 \pi \hat{\sigma}^2) - \frac{1}{2 \hat{\sigma}^2} (x -\hat \mu_k)^2 \\ &= \log(\hat{\pi}_k) -\frac{1}{2}\log(2 \pi \hat{\sigma}^2) - \frac{1}{2 \hat{\sigma}^2} \left(x^2 - 2 \hat \mu_k x + \mu_k^2 \right) \\ &= \log(\hat{\pi}_k) + x \frac{\hat{\mu}_k}{\hat{\sigma}^2} - \frac{\hat \mu_k^2}{2 \hat \sigma^2} -\frac{1}{2}\log(2 \pi \hat{\sigma}^2) - \frac{x^2}{2 \hat{\sigma}^2} \\ &= \log(\hat{\pi}_k) + x \frac{\hat{\mu}_k}{\hat{\sigma}^2} - \frac{\hat \mu_k^2}{2 \hat \sigma^2} + c \end{aligned} \]

Conclusion

Discriminant function discriminates based on \[ \log(\hat{\pi}_k) + x \frac{\hat{\mu}_k}{\hat{\sigma}^2} - \frac{\hat \mu_k^2}{2 \hat \sigma^2} \]

Question 2b

Show that if \(\hat{\mu}_1 < \hat{\mu}_2\), then every \(x \in \Re\) such that \(x < x^*\) is classified as class 1, and every \(x \in \Re\) such that \(x > x^*\) is classified as class 2, where \[ x^* = \frac{\hat{\mu}_1 + \hat{\mu}_2}{2} + \frac{\hat\sigma^2 \left[\log(\hat{\pi}_1) - \log(\hat{\pi}_2)\right]}{\hat{\mu}_2 - \hat{\mu}_1}. \]

Show

\[ \psi(x) = \mathbb{1}\{x < x^* \} + 1 \]

Know

• \(\psi(x) = \mathbb{1}\{\delta_1(x) > \delta_2(x)\} + 1\)
• \(\delta_k(x) = x \frac{\hat{\mu}_k}{\hat{\sigma}^2} - \frac{\hat{\mu}_k^2}{2 \hat{\sigma}^2} + \log(\hat{\pi}_k)\)

Question 2b

Plug in

\[ \begin{aligned} \delta_1(x) > \delta_2(x) &\iff x \frac{\hat{\mu}_1}{\hat{\sigma}^2} - \frac{\hat{\mu}_1^2}{2 \hat{\sigma}^2} + \log(\hat{\pi}_1) > x \frac{\hat{\mu}_2}{\hat{\sigma}^2} - \frac{\hat{\mu}_2^2}{2 \hat{\sigma}^2} + \log(\hat{\pi}_2) \\ &\iff \frac{\hat \mu_2^2 - \hat{\mu}_1^2}{2 \hat{\sigma}^2} + \log(\hat{\pi}_1) - \log(\hat{\pi}_2) > x \frac{\hat{\mu}_2-\hat{\mu}_1}{\hat{\sigma}^2} \\ &\iff \frac{1}{2 } \frac{\hat \mu_2^2 - \hat{\mu}_1^2}{\hat{\mu}_2-\hat{\mu}_1} +\hat{\sigma}^2\frac{[\log(\hat{\pi}_1) - \log(\hat{\pi}_2)]}{\hat{\mu}_2-\hat{\mu}_1} > x \\ &\iff \frac{1}{2 } \frac{(\hat \mu_2 - \hat{\mu}_1)(\hat \mu_2 + \hat{\mu}_1)}{\hat{\mu}_2-\hat{\mu}_1} +\hat{\sigma}^2\frac{[\log(\hat{\pi}_1) - \log(\hat{\pi}_2)]}{\hat{\mu}_2-\hat{\mu}_1} > x \\ &\iff \frac{1}{2 } \frac{\hat \mu_2 + \hat{\mu}_1}{\hat{\mu}_2-\hat{\mu}_1} +\hat{\sigma}^2\frac{[\log(\hat{\pi}_1) - \log(\hat{\pi}_2)]}{\hat{\mu}_2-\hat{\mu}_1} > x \\ &\iff x^* > x \end{aligned} \]

Question 2b

xs <- seq(from = -5, to = 5, length.out = 1000) 

mu_1 <- -0.5 
mu_2 <- 0.5

pi_1 <- 0.5
pi_2 <- 0.5 

sigma <- 1

x_star = (mu_1 + mu_2)/2 + sigma^2 * (log(pi_1) - log(pi_2)) / (mu_2 - mu_1) 

plot(xs, pi_1 * dnorm(xs, mean = mu_1, sd = sigma), type = "l", col = "green", lwd = 2, main = "Equal pis") 
lines(xs, pi_2 * dnorm(xs, mean = mu_2, sd = sigma), type = "l", col = "blue", lwd = 2)
abline(v = x_star, col = "grey", lty = 2) 
legend("topright", col = c("green", "blue", "grey"), lty = c(1,1,2), lwd = rep(2,3), legend = c("Class 1", "Class 2", "x_star")) 

pi_1 <- 0.6
pi_2 <- 0.4

x_star = (mu_1 + mu_2)/2 + sigma^2 * (log(pi_1) - log(pi_2)) / (mu_2 - mu_1) 


plot(xs, pi_1 * dnorm(xs, mean = mu_1, sd = sigma), type = "l", col = "green", lwd = 2, main = "Unequal pis") 
lines(xs, pi_2 * dnorm(xs, mean = mu_2, sd = sigma), type = "l", col = "blue", lwd = 2)
abline(v = x_star, col = "grey", lty = 2) 
legend("topright", col = c("green", "blue", "grey"), lty = c(1,1,2), lwd = rep(2,3), legend = c("Class 1", "Class 2", "x_star")) 

Learnings

• From definition to statement
• Under equal \(\hat\pi_k\), boundary is midway point between class means
• If \(\hat \pi_1 > \hat{\pi}_2\), second term positive, boundary shift to the right (favour class 1), and vice-versa

Question 3

Consider a three-class classification problem, with \(Y = i\) for class \(i \in \{1,2,3\}\).

Let \(X = (X_1, X_2)^\top\), with the following mean vectors:

\[ \mu_1 = (0, 0)^\top, \quad \mu_2 = (1, 1)^\top, \quad \mu_{3,1} = (0.5, 0.5)^\top, \quad \mu_{3,2} = (-0.5, 0.5)^\top \,, \]

and the covariance matrix is given by

\[ \Sigma = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\,. \]

• Suppose that \(X \mid Y = 1\) and \(X \mid Y = 2\) follow bivariate Gaussian distributions \(N(\mu_1, \Sigma)\) and \(N(\mu_2, \Sigma)\), respectively.

• Suppose \(X \mid Y = 3\) follows a mixture of two bivariate Gaussian distributions, with a 50% chance of being from \(N(\mu_{3,1}, \Sigma)\) and a 50% chance from \(N(\mu_{3,2}, \Sigma)\).

Additionally, assume

\[ \mathbb{P}[Y=1] = \tfrac{2}{5}, \quad \mathbb{P}[Y=2] = \tfrac{2}{5}, \quad \mathbb{P}[Y=3] = \tfrac{1}{5}. \]

Classify the point \(x_0 = (0.3, 0.3)^\top\) based on the posterior probabilities.

Recap

Posterior distribution

Bayes’ Theorem:

\[ \Pr(Y = k \mid X = x) = \frac{\pi_k f_k(x)}{\sum_{k = 1}^K \pi_k f_k(x)} \]

Classification:

\[ \underset{k}{\arg \max}\, \Pr(Y = k \mid X = x) = \underset{k}{\arg \max} \, \pi_k f_k(x) \]

Recap

Mixture distribution

A random variable \(X\) follows a mixture distribution with weights \(w_1, \ldots, w_C\) and mixture densities \(f_1,\ldots,f_C\) if:

\[ \Pr(X \leq t) = \sum_{c = 1}^k w_c F_c(t)\,, \]

where \(F_c(t) = \int_{-\infty}^t f_c(x) dx\) is the CDF associated with \(f_c\).

c1 <- dnorm(xs, mean = mu_1, sd = 1)
c2 <- ifelse(xs > 0, dexp(xs), 0)

plot(xs, c1, type = "l", lty = 2, col = "green", lwd = 2, main = "A mixture distribution", ylim = c(0,1)) 
lines(xs, c2, type = "l", lty = 2, col = "blue", lwd = 2)
lines(xs, pi_1 * c1 + pi_2 * c2, col = "orange", lwd = 2) 
legend("topright", col = c("green", "blue", "orange"), lty = c(2,2,1), lwd = rep(2,3), legend = c("Class 1", "Class 2", "Mixture")) 

Recap

Mixture distribution

Can think about a draw of \(X\) as follows:

1. Draw mixture class \(c\) from classes \(1,\ldots,C\) with probabilities \(w_1, \ldots, w_C\)
2. Draw \(X\) according to \(f_c\)

\[ \begin{aligned} \Pr (X \leq t) &= \Pr\left(\bigcup_{c = 1}^C \{\text{class } = c, X \leq t\}\right) \\ &= \sum_{c=1}^C \Pr( \text{class } = c, X \leq t) \\ &= \sum_{c=1}^C \Pr( \text{class } = c) \Pr( X \leq t \mid \text{class } = c) \\ &= \sum_{c=1}^C w_c F_c(t) \,. \end{aligned} \]

Recap

Mixture distribution

PDF of mixture distribution: \[ \begin{aligned} f(x) &= \frac{\partial}{\partial x} F(x) \\ &= \frac{\partial}{\partial x} \sum_{c = 1}^k w_c F_c(t) \\ &= \sum_{c = 1}^kw_c \frac{\partial}{\partial x} F_c(t) \\ &= \sum_{c = 1}^kw_c f_c(t) \\ \end{aligned} \]

Question 3

Show

How to classify \(x = (0.3,0.3)\) according to posterior probabilities

Know

• \(\Pr(Y = k \mid X = x) \propto \pi_k f_k\)
• \([\pi_1, \pi_2, \pi_3] = [0.4, 0.4, 0.2]\)
• The conditional pdfs of \(X\) conditional on label \(Y\) (two normal, one mixture)

Question 3

Plug-in

We have:

\[ \pi_1 f_1(x_0) = \frac{2}{5} \times \frac{1}{2\pi} \exp\!\left( -\frac{1}{2}(0.3 - 0,\, 0.3 - 0)(0.3 - 0,\, 0.3 - 0)^\top \right) = 0.0582, \]

\[ \pi_2 f_2(x_0) = \frac{2}{5} \times \frac{1}{2\pi} \exp\!\left( -\frac{1}{2}(0.3 - 1,\, 0.3 - 1)(0.3 - 1,\, 0.3 - 1)^\top \right) = 0.0390, \]

and

\[ \begin{aligned} \pi_3 f_3(x_0) &= \frac{1}{5} \times \frac{1}{2} \times \frac{1}{2\pi} \exp\!\left( -\frac{1}{2}(0.3 - 0.5,\, 0.3 - 0.5)(0.3 - 0.5,\, 0.3 - 0.5)^\top \right) \\ &\quad + \frac{1}{5} \times \frac{1}{2} \times \frac{1}{2\pi} \exp\!\left( -\frac{1}{2}(0.3 + 0.5,\, 0.3 + 0.5)(0.3 - 0.5,\, 0.3 - 0.5)^\top \right)\\ &= 0.0266. \end{aligned} \]

Since \(\pi_1 f_1(x_0)\) is the largest, we classify \(x_0\) as belonging to class 1.

Question 3

Learnings

• Generalise the discriminant approach \(\delta_k(x) = \log(\pi_k f_k)\) to unseen setting

Question 4

In this exercise, we consider binary classification. a. Consider a classifier that assigns each example to the majority label in the training dataset. What is the expected test accuracy of this classifier?

Show

Expected test accuracy of majority vote classifier

Know

• \(\psi(x) = \mathbb{1}\left\{\frac{1}{n}\sum_{i = 1}^n Y_i > \frac{1}{2} \right\} = \psi\)
• Accuracy of \(\psi\) at \(Y\): \(\mathbb{1}\{Y = \psi\}\)
• Expected accuracy: \(\mathrm{E}_{Y, \mathcal{D}}[\mathbb{1}\{Y = \psi\}]\) (\(\mathcal{D}\) is the training data with which we train our classifier)
• \(Y \perp \psi\)

Question 4a

Plug in

\[ \begin{aligned} \mathrm{E}_{Y, \mathcal{D}}[\mathbb{1}\{Y = \psi\}] &= \Pr\left(Y = \psi \right) \\ &= \Pr\left( \{Y = \psi = 1\} \cup \{Y = \psi = 0\} \right) \\ &= \sum_{i = 0}^1 \Pr \left( Y = \psi = i \right) \\ &= \sum_{i = 0}^1 \Pr \left( Y = i \right) \Pr(\psi = i) \\ &= \pi \pi_\psi + (1 - \pi) (1 - \pi_\psi) \,, \end{aligned} \] where \[ \pi = \Pr(Y = 1), \quad \pi_\psi = \Pr(\psi = 1) \]

Question 4a

What is \(\Pr(\psi = 1)\)?

To make further progress, we need to understand \(\psi\). Let us stick to our guns and write down what we know

\[ \Pr(\psi = 1) = \Pr\left(\mathbb{1}\left\{\frac{1}{n}\sum_{i = 1}^n Y_i > \frac{1}{2}\right\} = 1 \right) = \Pr \left( \frac{1}{n}\sum_{i = 1}^n Y_i > \frac{1}{2} \right) = \Pr \left( \sum_{i = 1}^n Y_i > \frac{n}{2} \right) \]

Now the random variable \(\sum_{i = 1}^n Y_i\) just counts the number of successes in \(n\) independent Bernoulli trials. Hence it follows a binomial distribution, and thus

\[ \begin{aligned} \Pr\left(\sum_{i = 1}^n Y_i > \frac{n}{2} \right) &= 1 - \Pr\left(\sum_{i = 1}^n Y_i \leq \frac{n}{2} \right) \\ &= 1 - \sum_{k = 1}^{\lfloor n / 2 \rfloor} {n\choose k} \pi^k (1- \pi)^{n-k} \end{aligned} \]

Question 4a

Hence we have found that

\[ \mathrm{E}_{Y, \mathcal{D}}[\mathbb{1}\{Y = \psi\}] = \pi \left(1 - \sum_{k = 1}^{\lfloor n / 2 \rfloor} {n\choose k} \pi^k (1- \pi)^{n-k} \right) + (1 - \pi) \left(\sum_{k = 1}^{\lfloor n / 2 \rfloor} {n\choose k} \pi^k (1- \pi)^{n-k} \right) \,. \]

Now although this is a precise characterisation of the expected prediction accuracy of our majority vote classifier, it looks rather convoluted and altogether a bit unsatisfying.

Intuitevly, we would expect that if we had enough votes, the majority vote should be \(1\) if the individual voting probabilities \(\pi\) are greater than \(1/2\) and \(0\) if they are less than \(1/2\).

Can we formalise this thought?

Question 4a

Intuitevly, we would expect that if we had enough votes, the majority vote should be \(1\) if the individual voting probabilities \(\pi\) are greater than \(1/2\) and \(0\) if they are less than \(1/2\).

n <- 10^3 

pi <- 0.7 

pi2 <- 0.3 

y1 <- rbinom(n, 1, pi) 
y2 <- rbinom(n, 1, pi2) 

av1 <- cumsum(y1) / 1:n 
av2 <- cumsum(y2) / 1:n 

plot(1:n, av1, type = "l", col = "blue", lwd = 2, main = "Average of responses", xlab = "Number of votes", ylab = "% of 1-votes") 
lines(1:n, av2, col = "orange", lwd = 2) 
abline(h = 0.5, col = "grey", lwd = 2, lty = 2) 
legend("bottomright", legend = c(expression(pi > 1/2), expression(pi < 1/2), "Decision boundary"), lty = c(1,1,2), lwd = rep(1,3), col = c("blue", "orange", "grey"))

Can we formalise this thought?

Recap

Law of large numbers

An average \(S_n = \sum_{i = 1}^n X_i / n\) of \(n\) i.i.d. random variables with expectation \(E[X] = \mu\) converges to \(\mu\) in probability as \(n \to \infty\): For all \(\epsilon, \delta\), there is \(N \in \mathbb{N}\), such that for all \(n > N\),

\[ \Pr(|S_n - \mu| > \epsilon ) < \delta \]

n <- 10^4 

pi <- 0.55 

eps <- 1e-2 
ub <- pi + eps
lb <- pi - eps 

set.seed(123)
y1 <- rbinom(n, 1, pi) 

av1 <- cumsum(y1) / 1:n 

plot(1:n, av1, type = "l", col = "blue", lwd = 2, main = "High probability bands", ylim = c(pi - 0.1, pi + .1), xlab = "Number of samples", ylab = expression(S[n])) 
abline(h = ub, col = "grey", lty = 2, lwd = 2) 
abline(h = lb, col = "grey", lty = 2, lwd = 2) 

polygon(
  x = c(par("usr")[1], par("usr")[2], par("usr")[2], par("usr")[1]),
  y = c(lb, lb, ub, ub),
  col = rgb(0, 0, 1, 0.2),
  border = NA
)

legend("bottomright", legend = c(expression(S[n]), expression(mu + epsilon), expression(mu - epsilon)), lty = c(1,2,2), lwd = rep(1,3), col = c("blue", "grey", "grey"))

Question 4a

\[ S_n = \frac{1}{n}\sum_{i = 1}^n Y_i \]

is an average of \(n\) i.i.d. Bernoulli random variables with success probability \(\pi\).

By the law of large numbers: \[ S_n \overset{p}{\to} \mathrm{E}[S_n] = \pi, \quad \text{as } n \to \infty \ \]

Question 4a

\(\pi > 1/2\)

n <- 10^3 

pi <- 0.6 

y1 <- rbinom(n, 1, pi) 

av1 <- cumsum(y1) / 1:n 

eps <- 1e-2 
ub <- pi + eps
lb <- pi - eps 



plot(1:n, av1, type = "l", col = "blue", lwd = 2, main = expression(pi > frac(1,2))) 
abline(h = ub, col = "grey", lty = 2, lwd = 2) 
abline(h = lb, col = "grey", lty = 2, lwd = 2) 
abline(h = 0.5, col = "red", lwd = 2, lty = 2) 
polygon(
  x = c(par("usr")[1], par("usr")[2], par("usr")[2], par("usr")[1]),
  y = c(lb, lb, ub, ub),
  col = rgb(0, 0, 1, 0.2),
  border = NA
)
polygon(
  x = c(par("usr")[1], par("usr")[2], par("usr")[2], par("usr")[1]),
  y = c(.5, .5, lb, lb),
  col = rgb(1, 0, 0, 0.2),  
  border = NA
)

legend("topright", legend = c(expression(S[n]), expression(pi + epsilon), expression(pi - epsilon), "Decision boundary"), lty = c(1,2,2,2), lwd = rep(1,4), col = c("blue", "grey", "grey", "red"))

• \(S_n \overset{p}{\to} \pi\)
• \(S_n - 1/2 \overset{p}{\to} \pi - 1/2 > 0\)
• \(\Pr(S_n - 1/2 > 0) \to 1\), as \(n \to \infty\)

Question 4a

\(\pi < 1/2\)

n <- 10^3 

pi <- 0.4 

set.seed(123)
y1 <- rbinom(n, 1, pi) 

av1 <- cumsum(y1) / 1:n 

eps <- 1e-2 
ub <- pi + eps
lb <- pi - eps 



plot(1:n, av1, type = "l", col = "blue", lwd = 2, main = expression(pi < frac(1,2))) 
abline(h = ub, col = "grey", lty = 2, lwd = 2) 
abline(h = lb, col = "grey", lty = 2, lwd = 2) 
abline(h = 0.5, col = "red", lwd = 2, lty = 2) 
polygon(
  x = c(par("usr")[1], par("usr")[2], par("usr")[2], par("usr")[1]),
  y = c(lb, lb, ub, ub),
  col = rgb(0, 0, 1, 0.2),
  border = NA
)
polygon(
  x = c(par("usr")[1], par("usr")[2], par("usr")[2], par("usr")[1]),
  y = c(ub, ub, 0.5, 0.5),
  col = rgb(1, 0, 0, 0.2),  
  border = NA
)

legend("bottomright", legend = c(expression(S[n]), expression(pi + epsilon), expression(pi - epsilon), "Decision boundary"), lty = c(1,2,2,2), lwd = rep(1,4), col = c("blue", "grey", "grey", "red"))

• \(\Pr(S_n - 1/2 > 0) \to 0\), as \(n \to \infty\)

Question 4a

\(\pi = 1/2\)

n <- 10^3 

pi <- 0.5

y1 <- rbinom(n, 1, pi) 

av1 <- cumsum(y1) / 1:n 

plot(1:n, av1, type = "l", col = "blue", lwd = 2, main = expression(pi==frac(1,2))) 
abline(h = 0.5, lwd = 2, lty = 2, col = "grey")

The process is \(S_n\) is symmetric. Each outcome \([y_1, \ldots, y_n]\) such that \(S_n > 1/2\) has a corresponding outcome \([1-2y_1, \ldots, 1-2y_n]\) with \(S_n < 1/2\) that has the same probability.

• \(\Pr(S_n > 1/2) = \Pr(S_n < 1 / 2)\)
• If \(n\) is odd: \(S_n = 1/2\) impossible: \(\Pr(S_n > 1/2) = 1/2\)
• If \(n\) is even: \(\Pr(S_n = 1/2) = (1/2)^n\). \(\Pr(S_n > 1/2) = (1 - (1/2)^n) / 2 \to 1/2\) as \(n \to \infty\)

Question 4a

\(\pi > 1/2\)

• \(\Pr(S_n > 1/2) = \Pr(\psi = 1) = \pi_\psi \to 1\) as \(n \to \infty\)
• \(\mathrm{E}_{Y, \mathcal{D}}[\mathbb{1}\{Y = \psi\}] = \pi \pi_\psi + (1 - \pi)(1 - \pi_\psi) \to \pi\)

\(\pi < 1/2\)

• \(\Pr(S_n > 1/2) = \Pr(\psi = 0) = 1 - \pi_\psi \to 1\) as \(n \to \infty\)
• \(\mathrm{E}_{Y, \mathcal{D}}[\mathbb{1}\{Y = \psi\}] = \pi \pi_\psi + (1 - \pi)(1 - \pi_\psi) \to 1 - \pi\)

\(\pi = 1/2\)

• \(\Pr(S_n > 1/2) = \pi_\psi = 1/2\)
• \(\mathrm{E}_{Y, \mathcal{D}}[\mathbb{1}\{Y = \psi\}] = \pi \pi_\psi + (1 - \pi)(1 - \pi_\psi) \to \pi = 1 - \pi\)

Question 4a

Overall: \[ \mathrm{E}_{Y, \mathcal{D}}[\mathbb{1}\{Y = \psi\}] \to \max\{\pi, 1-\pi\}, \quad \text{as } n \to \infty \]

Question 4a

Learnings

• Use definitions to characterise unknown classifier
• Use of rigour (definitions, LLN) and intuition (Majority class should be class with higher probability) to get a precise result
• Majority rate classifier classifies majority class correctly with probability going to one \((\pi \neq 1/2)\)
• Serves as a lower bound in accuracy for more sophisticated ML algorithms

Question 4b

For an email spam filter, it may be preferable to allow for some rate of spam to reach the inbox rather than risk sending an important legitimate email message to the spam folder. Assume we train a spam classifier where the positive class is spam and the negative class is non-spam.

Suppose we have three points on the ROC curve for our classifier: A: \((0.1, 0.7)\), B: \((0.2, 0.75)\) and C: \((0.3, 0.8)\). Which of these is preferable?

Recap

ROC

Binary, discriminant based classifiers are of the form \(\psi(x) = \mathbb{1}\{ s(x) < t \}\), where

• \(s(x)\) is a score that maps \(x\) to the real numbers
• \(t\) is treshold that translates scores into predictions

The scores induce an ordering on a test set \(X_1, \ldots, X_n\):

\[ s(X_{(1)}) < s(X_{(2)}) < \ldots < s(X_{(i)}) < t < s(X_{(i + 1)}) < \ldots < s(X_{(n)}) \]

where \(X_{(1)}\) is the feature that produces the lowest score and so on.

• If \(t < s(X_{(1)})\), then all points are labelled zero and both the TPR and FPR are zero
• For \(s(X_{(i+1)}) \geq t > s(X_{(i)})\), \(X_{(1)}, \ldots, X_{(i)}\) are labelled positive, the rest negative
• ROC plots \(TPR(t), FPR(t)\) for all values of \(t\)

The ROC curve gives an indication of how well a prediction method discriminates among classes, regardless of the threshold

Recap

Confusion matrix for binary responses

\[ \text{TPR} = \frac{TP}{TP + FN} \]

\[ \text{FPR} = \frac{FP}{FP + TN} \]

Recap

ROC curve

# Load necessary libraries
library(dplyr)
library(knitr)
library(kableExtra) 

# Set seed for reproducibility
set.seed(123)

# Define the number of observations
n <- 70

# Generate random data for each covariate

# Age: Assuming cat age follows an absolute normal distribution
age <- abs(rnorm(n, sd = 4))

# Color: Ginger, Black, Brown
color <- c("Ginger", "Black", "Brown")[rbinom(n, 2, age/max(age)) + 1]

# Multiple Colors: binary with conditional probabilities
multiple_colors <- sapply(color, function(col) {
  if (col == "Ginger") {
    return(c("No", "Yes")[rbinom(1, 1, 0.6) + 1]) 
  } else if (col == "Black") {
    return(c("No", "Yes")[rbinom(1, 1, 0.05) + 1]) 
  } else {
    return(c("No", "Yes")[rbinom(1, 1, 0.25) + 1]) 
  }
})

# Generate previous_owners correlated with age
previous_owners <- sapply(age, function(a) {
  prob <- c(0.05, 0.7, 0.2, 0.05) * (a / max(age))
  prob <- prob / sum(prob) 
  return(sample(0:3, 1, prob = prob))
})
# Neutered: Binary with p = 0.9
neutered <- c("No","Yes")[rbinom(n, 1, 0.9) + 1]

# Gender: Binary with p = 0.5
gender <- c("Male", "Female")[rbinom(n, 1, ifelse(multiple_colors == "Yes", 0.1, 0.5)) + 1]

# Playfulness: range between 1 and 10, negatively correlated with age and col == black, positively with neutered and gender 
playfulness <- -0.2 * age + rnorm(n, sd = 2) + 0.25 * (neutered == "Yes") - 0.1 * (2 * (gender == "Male") - 1) - 5 * (color == "Black")
playfulness <- 10 * ((playfulness - min(playfulness)) / (max(playfulness) - min(playfulness)))

# Create a data frame
cat_adoption_data <- data.frame(
  Age = round(age, 1),
  MainColor = as.factor(color),
  MultipleColors = as.factor(multiple_colors),
  PreviousOwners = previous_owners,
  Neutered = as.factor(neutered),
  Gender = as.factor(gender),
  Playfulness = round(playfulness, 1)
)

# Expand categorical variables using model.matrix
model_data <- model.matrix(~ Age + MainColor + MultipleColors + PreviousOwners + Neutered + Gender + Playfulness, data = cat_adoption_data)

# Define a true effect vector (for simplicity, assume some arbitrary effects)
true_effect_base <- c(Age = -0.5, MainColorBrown = .3, MainColorGinger = -.5, MultipleColorsYes = 0.25, PreviousOwners = 0, NeuteredYes = 0, GenderMale = 0, Playfulness = 0.25)

true_effects <- c(-c(Intercept = mean(model_data[,-1] %*% true_effect_base)), true_effect_base) 

# Simulate the probability of adoption using the logistic function
prob_adoption <- plogis(model_data %*% true_effects)

# Simulate the binary adoption outcome
adopted <- c("No", "Yes")[rbinom(n, 1, prob_adoption) + 1]

# Add the adoption outcome to the data frame
cat_adoption_data <- cbind(cat_adoption_data, Adopted = adopted)

Question 4b

Show

• Which point on ROC curve is most preferrable for a spam filter

Know

• ROC points record \(TPR(t), FPR(t)\) for a particular threshold \(t\)
• TP: Spam emails classified as spam
• FP: Normal emails classified as spam
• Prefer spam emails classified as normal (FN), over normal emails classified as spam (FP)
• FNR = 1 - TPR

Question 4b

Look for point with low FPR, at expense of lower TPR

x0 <- c(0,0)
x1 <- c(0.1, 0.7) 
x2 <- c(.2, 0.75) 
x3 <- c(.3, .8)
x4 <- c(1,1)

xs <- rbind(x0, x1, x2, x3, x4) 

plot(xs, type = "b", col = "blue", lwd = 2, xlab = "FPR", ylab = "TPR", main = "ROC") 
abline(a = 0, b = 1, col = "grey", lwd = 2, lty = 2) 

ROC_points <- rbind(x1,x2,x3) 
point_labels <- c("A", "B", "C") 
text(ROC_points, point_labels, pos=3, cex = 1.5)

• A has lowest FPR, with comparable TPR, so it is most preferable