hw3.Rmd

---
title: "Homework 3"
output: html_document
username: aamster3
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(ggplot2)
```

# Data Preprocessing
```{r}
train = t(read.csv('mnist/mnist_train.csv', header = F))
test = t(read.csv('mnist/mnist_test.csv', header = F))

train_0_1 = train[train[, 785] == 0 | train[, 785] == 1, ]
train_3_5 = train[train[, 785] == 3 | train[, 785] == 5, ]

test_0_1 = test[test[, 785] == 0 | test[, 785] == 1, ]
test_3_5 = test[test[, 785] == 3 | test[, 785] == 5, ]

print(dim(train_0_1))
print(dim(train_3_5))
print(dim(test_0_1))
print(dim(test_3_5))

train_labels_0_1 = as.vector(train_0_1[, 785])
test_labels_0_1 = as.vector(test_0_1[, 785])
train_labels_3_5 = as.vector(train_3_5[, 785])
test_labels_3_5 = as.vector(test_3_5[, 785])

train_0_1 = train_0_1[, 1:784]
test_0_1 = test_0_1[, 1:784]
train_3_5 = train_3_5[, 1:784]
test_3_5 = test_3_5[, 1:784]
```
```{r}
rotate = function(x) t(apply(x, 2, rev))

index0 = which(train_labels_0_1 == 0)[1]
index1 = which(train_labels_0_1 == 1)[1]
index3 = which(train_labels_3_5 == 3)[1]
index5 = which(train_labels_3_5 == 5)[1]
```
```{r}
print(paste('label: ', as.character(train_labels_0_1[index0])))
image0 = matrix(train_0_1[index0, ], nrow = 28, ncol = 28)
image(rotate(image0), col=gray(0:255/255))
```
```{r}
print(paste('label: ', train_labels_0_1[index1]))
image1 = matrix(train_0_1[index1, ], nrow = 28, ncol = 28)
image(rotate(image1), col=gray(0:255/255))
```
```{r}
print(paste('label: ', train_labels_3_5[index3]))
image3 = matrix(train_3_5[index3, ], nrow = 28, ncol = 28)
image(rotate(image3), col=gray(0:255/255))
```
```{r}
print(paste('label: ', train_labels_3_5[index5]))
image5 = matrix(train_3_5[index5, ], nrow = 28, ncol = 28)
image(rotate(image5), col=gray(0:255/255))
```
# Implementation
```{r}
DEBUG = F

# Set labels in train_0_1 to be 1 if 0 else 0
# Set labels in train_3_5 to be 1 if 3 else 0
# do the same for test data

train_labels_0_1 = ifelse(train_labels_0_1 == 0,  1, 0)
train_labels_3_5 = ifelse(train_labels_3_5 == 3,  1, 0)
test_labels_0_1 = ifelse(test_labels_0_1 == 0,  1, 0)
test_labels_3_5 = ifelse(test_labels_3_5 == 3,  1, 0)

# add bias term to training and testing data

train_0_1 = cbind(rep(1, nrow(train_0_1)), train_0_1)
train_3_5 = cbind(rep(1, nrow(train_3_5)), train_3_5)
test_0_1 = cbind(rep(1, nrow(test_0_1)), test_0_1)
test_3_5 = cbind(rep(1, nrow(test_3_5)), test_3_5)

g = function(z) {
  #print(paste('g', 1 / (1 + exp(-z))))
  return(1 / (1 + exp(-z)))
}

h = function(theta, x) {
  #print(paste('z', t(theta) %*% x))
  #print(paste('g(z)', g(c(t(theta) %*% x))))
  return(g(c(t(theta) %*% x)))
}

cost = function(yhat, y) {
  #print(paste('y', y))
  #print(paste('yhat', yhat))
  return(mean(-(y * log(yhat) + (1 - y) * log(1 - yhat))))
}

train = function(data, labels, alpha = 1e-3, epsilon = 1e-4, max_iterations = 20) {
  theta = rnorm(ncol(data)) 
  
  converged = F
  
  
  for (iter in 1:max_iterations) {
    sample_indices = sample(1:nrow(data))
    data = data[sample_indices, ]
    labels = labels[sample_indices]
    
    m = nrow(data)
    for (i in 1:m) {
      x = data[i, ]
      y = labels[i]
      gradient = (h(theta, x) - y) * x
      update = alpha * gradient
      theta = theta - update
    }
    
    gradient_norm = norm(gradient, "2")
    
    if (gradient_norm < epsilon) { 
      converged = T
      break
    }
    
    if (DEBUG) {
     print(paste('accuracy', accuracy(predict(theta, data), labels)))
      
      if (!converged){
        print(paste('max iter reached', 'gradient norm', gradient_norm))
      }
    }
  }
  
  return(theta)
}

predict = function(theta, data) {
  pred = apply(data, 1, function(x) h(theta, x))
  return(ifelse(pred > 0.5, 1, 0))
}

accuracy = function(pred, labels) {
  return(sum(pred == labels) / length(labels))
}

print('model 0_1')
model_0_1 = train(train_0_1, train_labels_0_1)

print('model 3_5')
model_3_5 = train(train_3_5, train_labels_3_5)
```
I initialize the coefficients randomly by sampling from the normal distribution.
The convergence criteria is when the l2 norm (▽theta) < epsilon
I took the transpose of the data in order to have the observations along the rows and the labels as a column vector. I added a bias term by inserting a column of 1s in the training data.

```{r}
# model 0_1

# 2 correct
pred = predict(model_0_1, train_0_1)
correct = which(pred == train_labels_0_1)
correct_sample_idx = correct[sample(1:length(correct), 2)]

# 2 incorrect
incorrect = which(pred != train_labels_0_1)
incorrect_sample_idx = incorrect[sample(1:length(incorrect), 2)]

# Drop the bias term for plotting
no_bias_train_0_1 = train_0_1[, 2:ncol(train_0_1)]
```
## 2 correct followed by 2 incorrect from train_0_1
```{r}

print(paste('pred', pred[correct_sample_idx[1]]))
print(paste('label', train_labels_0_1[correct_sample_idx[1]]))
image = matrix(no_bias_train_0_1[correct_sample_idx[1], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```

```{r}
print(paste('pred', pred[correct_sample_idx[2]]))
print(paste('label', train_labels_0_1[correct_sample_idx[2]]))
image = matrix(no_bias_train_0_1[correct_sample_idx[2], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```
```{r}
print(paste('pred', pred[incorrect_sample_idx[1]]))
print(paste('label', train_labels_0_1[incorrect_sample_idx[1]]))
image = matrix(no_bias_train_0_1[incorrect_sample_idx[1], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```

```{r}
print(paste('pred', pred[incorrect_sample_idx[2]]))
print(paste('label', train_labels_0_1[incorrect_sample_idx[2]]))
image = matrix(no_bias_train_0_1[incorrect_sample_idx[2], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```
## 2 correct followed by 2 incorrect from train_3_5
```{r}
pred = predict(model_3_5, train_3_5)
correct = which(pred == train_labels_3_5)
correct_sample_idx = correct[sample(1:length(correct), 2)]

incorrect = which(pred != train_labels_3_5)
incorrect_sample_idx = incorrect[sample(1:length(incorrect), 2)]

# Drop the bias term for plotting
no_bias_train_3_5 = train_3_5[, 2:ncol(train_3_5)]
```
```{r}

print(paste('pred', pred[correct_sample_idx[1]]))
print(paste('label', train_labels_3_5[correct_sample_idx[1]]))
image = matrix(no_bias_train_3_5[correct_sample_idx[1], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```

```{r}
print(paste('pred', pred[correct_sample_idx[2]]))
print(paste('label', train_labels_3_5[correct_sample_idx[2]]))
image = matrix(no_bias_train_3_5[correct_sample_idx[2], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```
```{r}
print(paste('pred', pred[incorrect_sample_idx[1]]))
print(paste('label', train_labels_3_5[incorrect_sample_idx[1]]))
image = matrix(no_bias_train_3_5[incorrect_sample_idx[1], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```

```{r}
print(paste('pred', pred[incorrect_sample_idx[2]]))
print(paste('label', train_labels_3_5[incorrect_sample_idx[2]]))
image = matrix(no_bias_train_3_5[incorrect_sample_idx[2], ], nrow = 28, ncol = 28)
image(rotate(image), col=gray(0:255/255))
```
# Modeling
```{r}
accuracy = function(labels, labels_pred) {
  return(sum(labels_pred == labels) / length(labels))
}

model = function(train_data, train_labels, test_data, test_labels, alpha) {
  theta = train(train_data, train_labels, alpha)
  train_accuracy = accuracy(train_labels, predict(theta, train_data))
  test_accuracy = accuracy(test_labels, predict(theta, test_data))
  
  return(list(theta = theta, train_acc = train_accuracy, test_acc = test_accuracy))
}

results = function(train_data, train_labels, test_data, test_labels) {
  ALPHA = c(.001, .01, .02, .04, .08, .1)
  NUM_ITERATIONS = 10
  
  results = data.frame(alpha = numeric(0), mean = numeric(0), se = numeric(0), dataset = character(0))
  
  for (alpha in ALPHA) {
    train_acc_sum = 0
    train_acc_sum_2 = 0
    test_acc_sum = 0
    test_acc_sum_2 = 0
    
    for (iter in 1:NUM_ITERATIONS) {
      model = model(train_data, train_labels, test_data, test_labels, alpha)
      train_acc = model[['train_acc']]
      test_acc = model[['test_acc']]
      
      train_acc_sum = train_acc_sum + train_acc
      train_acc_sum_2 = train_acc_sum_2 + train_acc^2
      test_acc_sum = test_acc_sum + test_acc
      test_acc_sum_2 = test_acc_sum_2 + test_acc^2
    }
    
    train_acc_mean = train_acc_sum / NUM_ITERATIONS
    train_acc_std = sqrt(train_acc_sum_2 / NUM_ITERATIONS - train_acc_mean^2) 
    test_acc_mean = test_acc_sum / NUM_ITERATIONS
    test_acc_std = sqrt(test_acc_sum_2 / NUM_ITERATIONS - test_acc_mean^2)
    
    train_acc_standard_error = train_acc_std / sqrt(NUM_ITERATIONS)
    test_acc_standard_error = test_acc_std / sqrt(NUM_ITERATIONS)
    
    results = rbind(results, data.frame(alpha = alpha, mean = train_acc_mean, se = test_acc_standard_error, dataset = 'train'))
    results = rbind(results, data.frame(alpha = alpha, mean = test_acc_mean, se = test_acc_standard_error, dataset = 'test'))
  }
  
  return(results)
}

results_0_1 = results(train_0_1, train_labels_0_1, test_0_1, test_labels_0_1)
results_3_5 = results(train_3_5, train_labels_3_5, test_3_5, test_labels_3_5)
```
```{r}
ggplot(results_0_1, aes(x = alpha, y = mean, colour=dataset, group=dataset)) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), colour="black") +
  geom_line() +
  geom_point(size=3, shape=21, fill="white") +
  ylab('mean accuracy') +
  ggtitle('0/1 Dataset')
```


```{r}
ggplot(results_3_5, aes(x = alpha, y = mean, colour=dataset, group=dataset)) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), colour="black") +
  geom_line() +
  geom_point(size=3, shape=21, fill="white") +
  ylab('mean accuracy') +
  ggtitle('3/5 dataset')
```

The best test accuracy for the 0/1 dataset is ~.997 while the best test accuracy for the 3/5 dataset is ~.955. This slight difference in accuracy I think is because a 3 and 5 look more similar than a 0 and 1 so the model has a more difficult time differentiating between the two digits.

In order to implemenent a multiclass classifier, we can use the one-vs-rest approach. Essentially a separate model is trained where one class becomes 1 and the rest become 0. A model is trained for each class we are trying to predict. Then we can take each trained model and choose the maximum value from all the models and use that to predict the class.

# Learning Curves
```{r}
results = function(train_data, train_labels, test_data, test_labels) {
  TRAINING_EXAMPLE_FRACTION = seq(from = .01, to = 1, length.out = 10)
  NUM_ITERATIONS = 10
  ALPHA = 0.04
  
  results = data.frame(training_fraction = numeric(0), mean = numeric(0), se = numeric(0), dataset = character(0))
  
  for (training_examples_fraction in TRAINING_EXAMPLE_FRACTION) {
    sample_idx = sample(1:nrow(train_data), size = nrow(train_data) * training_examples_fraction)
    train_subset = train_data[sample_idx, ]
    train_subset_labels = train_labels[sample_idx]
    
    train_acc_sum = 0
    train_acc_sum_2 = 0
    test_acc_sum = 0
    test_acc_sum_2 = 0
    
    for (iter in 1:NUM_ITERATIONS) {
      model = model(train_subset, train_subset_labels, test_data, test_labels, ALPHA)
      train_acc = model[['train_acc']]
      test_acc = model[['test_acc']]
      
      train_acc_sum = train_acc_sum + train_acc
      train_acc_sum_2 = train_acc_sum_2 + train_acc^2
      test_acc_sum = test_acc_sum + test_acc
      test_acc_sum_2 = test_acc_sum_2 + test_acc^2
    }
    
    train_acc_mean = train_acc_sum / NUM_ITERATIONS
    train_acc_std = sqrt(train_acc_sum_2 / NUM_ITERATIONS - train_acc_mean^2) 
    test_acc_mean = test_acc_sum / NUM_ITERATIONS
    test_acc_std = sqrt(test_acc_sum_2 / NUM_ITERATIONS - test_acc_mean^2)
    
    train_acc_standard_error = train_acc_std / sqrt(NUM_ITERATIONS)
    test_acc_standard_error = test_acc_std / sqrt(NUM_ITERATIONS)
    
    results = rbind(results, data.frame(training_fraction = training_examples_fraction, mean = train_acc_mean, se = test_acc_standard_error, dataset = 'train'))
    results = rbind(results, data.frame(training_fraction = training_examples_fraction, mean = test_acc_mean, se = test_acc_standard_error, dataset = 'test'))
  }
  
  return(results)
}

results_0_1 = results(train_0_1, train_labels_0_1, test_0_1, test_labels_0_1)
results_3_5 = results(train_3_5, train_labels_3_5, test_3_5, test_labels_3_5)
```

```{r}
ggplot(results_0_1, aes(x = training_fraction, y = mean, colour=dataset, group=dataset)) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), colour="black") +
  geom_line() +
  geom_point(size=3, shape=21, fill="white") +
  ylab('mean accuracy') +
  ggtitle('0/1 Dataset')
```

```{r}
ggplot(results_3_5, aes(x = training_fraction, y = mean, colour=dataset, group=dataset)) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), colour="black") +
  geom_line() +
  geom_point(size=3, shape=21, fill="white") +
  ylab('mean accuracy') +
  ggtitle('3/5 Dataset')
```

In the above two plots we can see that the training and testing accuracies increase as we increase the training fraction from .01 to 1. We can see that the train accuracy and test accuracy are quite close to each other, a sign of high bias. The test accuracy essentially plateaus, meaning that this is about the best we can do using the given model. We should introduce parameters into the model to fit the data more closely or use a different machine learning model, then we may be able to increase test accuracy a little bit.