Mini Project 6

Watchanan Chantapakul (wcgzm)

---

This miniproject has only one question, which is a bit longer than the questions in previous assignments, however it also builds on the previous assignments. So, at the end of the day, you should be able to "borrow" a lot of your own code from before and finish this assignment quite easily. In this experiment, you are to compare various different classification approaches:

1. using complete knowledge of the statistics of the data and computing optimum discriminant functions for the classification (Chapter 2);
2. assuming that you know only the model of the distributions, but not their parameters (Chapter 3);
3. by first reducing the dimensionality of the data set and then classifying it (Chapter 3);
4. assuming that you do not know the underlying distributions and employing Parzen Window (Chapter 4);
5. assuming that you do not know the underlying distributions and employing k-NN (Chapter 4);
6. assuming that you do not know anything, but the class labels (Chapter 5).

---

%matplotlib inline

import os

import numpy as np
# save np.load
np_load_old = np.load
# modify the default parameters of np.load
np.load = lambda *a,**k: np_load_old(*a, allow_pickle=True, **k)

import scipy
import scipy.io

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.lines import Line2D
from mpl_toolkits.mplot3d import Axes3D

from IPython.display import display, Math

np.random.seed(7720)

def show(text, ans, precision=4):
    if type(ans) == np.ndarray:
        t = r'\begin{bmatrix} '
        for i in ans:
#             print(i, type(i))
#             print(r' \\ '.join(i))
            if type(i) != np.ndarray:
                t += f'{i:.{precision}f}' + r' \\ '
            else:
                a_str = np.array2string(i, precision=precision, separator=r' & ')
                t += a_str[1:-1]
                t += r' \\ '
        t += r'\end{bmatrix}'
        display(Math(f'{text} = {t}'))
    else:
        display(Math(f'{text} = {ans:.{precision}f}'))

def show_percent(text, ans, precision=2):
    display(Math(f'{text} = {ans:.{precision}f}\%'))

Generate dataset and save it to 'data/mp6.npy' file

Initially, you must create four data sets according to the information below, but next week, you will be given the actual datasets with which you will write your final reports. The datasets that you will create are 'trivial', but they will be useful for debugging your programs. The ones that I will provide will be more challenging.

While using your datasets, the four testing/training data sets described here MUST be kept the same through out all parts below. That is, you should create your data points ONCE, save them, and use them for all parts and subparts below. Do NOT create different data for each question or you may end up with very strange results!!

The four data sets will be referred to as: 1) Training Data I, with 50 samples in each class; 2) Training Data II, with 500 samples in each class; 3) Testing Data I, also with 500 samples in each class; and finally 4) Testing Data II, with 10000 samples in each class¹

All data sets must consist of 5-d points divided in 3-classes with the underlying normal distributions $p(\vec{x} | \omega_i) = \mathsf{N}(\vec{\mu}_i, \Sigma_i)$, where:

$$ \begin{align*} \vec{\mu}_1 &= \begin{bmatrix}2 & 3 & 1 & 5.5 & 8.7\end{bmatrix}^t\\ \vec{\mu}_2 &= \begin{bmatrix}-4.5 & 6 & -1 & 3 & 10\end{bmatrix}^t\\ \vec{\mu}_3 &= \begin{bmatrix}1.2, -2.3, 1.5, -0.5, 2.7\end{bmatrix}^t\\ \end{align*} $$$$ \vec{\Sigma}_1 = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 0.5 & 0 & 0 & 0 \\ 0 & 0 & 2.5 & 0 & 0 \\ 0 & 0 & 0 & 0.7 & 0 \\ 0 & 0 & 0 & 0 & 3.5 \\ \end{bmatrix} \vec{\Sigma}_2 = \begin{bmatrix} 2 & 0 & 1 & 0.5 & 0 \\ 0 & 3.5 & 0 & 0 & 0.6 \\ 1 & 0 & 4.5 & 1.2 & 0 \\ 0.5 & 0 & 1.2 & 1.6 & 0 \\ 0 & 0.6 & 0 & 0 & 2.5 \\ \end{bmatrix} \vec{\Sigma}_3 = \begin{bmatrix} 4.2 & 0 & 1.3 & 2.5 & 1.4 \\ 0 & 5 & 0 & 0 & 3.6 \\ 1.3 & 0 & 4.5 & 4.2 & 0 \\ 2.5 & 0 & 4.2 & 5.6 & 0 \\ 1.4 & 3.6 & 0 & 0 & 7.5 \\ \end{bmatrix} $$

You will assume that all states of nature are equally probable.

---

¹ The sizes of the data sets above are not typical to real-life classification problems. The choices above were made solely to create interesting situations for discussion in your report.

---

* Comment out the code for generating the dataset (run once for the first time)

# C = 3
# data = {}
# data['train1'] = [None] * C
# data['train2'] = [None] * C
# data['test1'] = [None] * C
# data['test2'] = [None] * C

# mu = [
#     np.array([2, 3, 1, 5.5, 8.7]),
#     np.array([-4.5, 6, -1, 3, 10]),
#     np.array([1.2, -2.3, 1.5, -0.5, 2.7]),
# ]
      
# cov = [
#     np.identity(5) * np.array([1, 0.5, 2.5, 0.7, 3.5]),
#     np.array([
#         [2, 0, 1, 0.5, 0],
#         [0, 3.5, 0, 0, 0.6],
#         [1, 0, 4.5, 1.2, 0],
#         [0.5, 0, 1.2, 1.6, 0],
#         [0, 0.6, 0, 0, 2.5]
#     ]),
#     np.array([
#         [4.2, 0, 1.3, 2.5, 1.4],
#         [0, 5, 0, 0, 3.6],
#         [1.3, 0, 4.5, 4.2, 0],
#         [2.5, 0, 4.2, 5.6, 0],
#         [1.4, 3.6, 0, 0, 7.5]
#     ]),
# ]    

# for c in range(C):
#     data['train1'][c] = np.random.multivariate_normal(mu[c], cov[c], 50)
#     data['train2'][c] = np.random.multivariate_normal(mu[c], cov[c], 500)
#     data['test1'][c] = np.random.multivariate_normal(mu[c], cov[c], 500)
#     data['test2'][c] = np.random.multivariate_normal(mu[c], cov[c], 10000)
    
# data['mu'] = mu
# data['cov'] = cov
# np.save('data/mp6.npy', data)

Load dataset from 'data/mp6.npy' file

data = np.load('data/mp6.npy').item()
C = len(data['mu'])

Cross-check the defined mean vectors and covariance matrices with the question.

show(r'\vec{\mu}_1^{\mathsf{T}}', data['mu'][0][:, None].T, precision=1)
show(r'\vec{\mu}_2^{\mathsf{T}}', data['mu'][1][:, None].T, precision=1)
show(r'\vec{\mu}_3^{\mathsf{T}}', data['mu'][2][:, None].T, precision=1)

$\displaystyle \vec{\mu}_1^{\mathsf{T}} = \begin{bmatrix} 2. & 3. & 1. & 5.5 & 8.7 \\ \end{bmatrix}$

$\displaystyle \vec{\mu}_2^{\mathsf{T}} = \begin{bmatrix} -4.5 & 6. & -1. & 3. & 10. \\ \end{bmatrix}$

$\displaystyle \vec{\mu}_3^{\mathsf{T}} = \begin{bmatrix} 1.2 & -2.3 & 1.5 & -0.5 & 2.7 \\ \end{bmatrix}$

show(r'\Sigma_1', data['cov'][0], precision=1)
show(r'\Sigma_2', data['cov'][1], precision=1)
show(r'\Sigma_3', data['cov'][2], precision=1)

$\displaystyle \Sigma_1 = \begin{bmatrix} 1. & 0. & 0. & 0. & 0. \\ 0. & 0.5 & 0. & 0. & 0. \\ 0. & 0. & 2.5 & 0. & 0. \\ 0. & 0. & 0. & 0.7 & 0. \\ 0. & 0. & 0. & 0. & 3.5 \\ \end{bmatrix}$

$\displaystyle \Sigma_2 = \begin{bmatrix} 2. & 0. & 1. & 0.5 & 0. \\ 0. & 3.5 & 0. & 0. & 0.6 \\ 1. & 0. & 4.5 & 1.2 & 0. \\ 0.5 & 0. & 1.2 & 1.6 & 0. \\ 0. & 0.6 & 0. & 0. & 2.5 \\ \end{bmatrix}$

$\displaystyle \Sigma_3 = \begin{bmatrix} 4.2 & 0. & 1.3 & 2.5 & 1.4 \\ 0. & 5. & 0. & 0. & 3.6 \\ 1.3 & 0. & 4.5 & 4.2 & 0. \\ 2.5 & 0. & 4.2 & 5.6 & 0. \\ 1.4 & 3.6 & 0. & 0. & 7.5 \\ \end{bmatrix}$

Confusion Matrix

All confusion matrices in this miniproject has the following structure:

	(predicted) 1	(predicted) 2	(predicted) 3
(actual) 1	-	-	-
(actual) 2	-	-	-
(actual) 3	-	-	-

Table 1: Confusion matrix format with the actual classes in rows and the predicted classes in columns.

Part I

Here, you will first (in a and b) use complete knowledge about the data. Then (in c and d), you will "forget" that you know the means and covariances, and use ML to estimate them.

a) Classify the Testing Data I, using the given statistics above and the Bayes decision rule. Compute the confusion matrix.

We have a priori knowledge about the distributions that they are Gaussian distributions. So, we can define the discriminant function for the normal density with the following generic form: $$ g_i(\vec{x}) = -\frac{1}{2}(\vec{x}-\vec{\mu}_i)^{\mathsf{T}}\mathbf{\Sigma}_i^{-1}(\vec{x}-\vec{\mu}_i)-\frac{d}{2}\ln (2\pi) - \frac{1}{2} \ln |\mathbf{\Sigma}_i| + \ln P(\omega_i) $$ for any given $d$-dimensional data, mean $\vec{\mu}_i$, covariance matrix $\mathbf{\Sigma}_i$ and prior probabilities $P(\omega_i)$ of class $i$.

Given the equal number of samples in every class, thus the priors are equal. We then can discard the term $\ln P(\omega_i)$. Also, since $d$ is a constant, we then also get rid of the constant term $- \frac{d}{2} \ln (2\pi)$.

We classify a sample $\vec{x}$, based on the Bayes decision rule, to be class $\omega_i$ if $g_i(\vec{x}) > g_j(\vec{x})$, $\forall j \neq i$. This can also be written in the form of argmax as follows: $$ \hat{\omega} = \underset{i}{\mathrm{argmax}} g_i(\vec{x}) $$

def squared_mahalanobis_distance(x, y, cov):
    a = np.array(x) - np.array(y)
    r2 = a.T @ np.linalg.inv(cov) @ a
    return r2

def discriminant_function(x, mean, cov):
    d = len(x)
    A = -0.5 * squared_mahalanobis_distance(x, mean, cov)
#     B = - ((d/2) * np.log(2 * np.pi))
    C = - (0.5 * np.log(np.linalg.det(cov)))
#     D = np.log(prior)
    return A + C

a_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test1'][c]:
        distances = np.zeros(C)
        for i in range(C):
            distances[i] = discriminant_function(x, data['mu'][i], data['cov'][i])
        pred_class = np.argmax(distances)
        a_confusion_matrix[c][pred_class] += 1

show('\mathrm{[Part\ I]\ a)\ \ Confusion\ matrix}', a_confusion_matrix)

$\displaystyle \mathrm{[Part\ I]\ a)\ \ Confusion\ matrix} = \begin{bmatrix} 499 & 1 & 0 \\ 1 & 499 & 0 \\ 0 & 0 & 500 \\ \end{bmatrix}$

def accuracy(cm):
    return cm.trace() / cm.sum(axis=None) * 100

show_percent(r'\mathrm{[Part\ I]\ a)\ accuracy}', accuracy(a_confusion_matrix))
show_percent(r'\mathrm{[Part\ I]\ a)\ error}', 100 - accuracy(a_confusion_matrix))

$\displaystyle \mathrm{[Part\ I]\ a)\ accuracy} = 99.87\%$

$\displaystyle \mathrm{[Part\ I]\ a)\ error} = 0.13\%$

b) Repeat part a) for the Testing Data II.

b_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        distances = np.zeros(C)
        for i in range(C):
            distances[i] = discriminant_function(x, data['mu'][i], data['cov'][i])
        pred_class = np.argmax(distances)
        b_confusion_matrix[c][pred_class] += 1

show('\mathrm{[Part\ I]\ b)\ \ Confusion\ matrix}', b_confusion_matrix)

$\displaystyle \mathrm{[Part\ I]\ b)\ \ Confusion\ matrix} = \begin{bmatrix} 9983 & 7 & 10 \\ 8 & 9992 & 0 \\ 9 & 0 & 9991 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ I]\ b)\ accuracy}', accuracy(b_confusion_matrix))
show_percent(r'\mathrm{[Part\ I]\ b)\ error}', 100 - accuracy(b_confusion_matrix))

$\displaystyle \mathrm{[Part\ I]\ b)\ accuracy} = 99.89\%$

$\displaystyle \mathrm{[Part\ I]\ b)\ error} = 0.11\%$

c) Compute the ML estimates $(\hat{\vec{\mu}}_i$ and $\hat{\Sigma}_i$) for each class using the Training Data I and classify the Testing Data II using Bayes decision rule. Compute the confusion matrix.

According to the maximum likelihood estimation (MLE) for a Gaussian distribution, the estimated mean vector of class $\omega_i$ can be computed using $$ \vec{\mu}_i = \frac{1}{n_i} \sum_{k=1}^{n_i} \vec{x}_k $$ where $n$ is the number of training samples in class $\omega_i$.

The estimated unbiased covariance matrix of class $\omega_i$ is given by: $$ \mathbf{\Sigma}_i = \frac{1}{n_i-1} \sum_{k=1}^{n_i} (\vec{x}_k - \vec{\mu}_i) (\vec{x}_k - \vec{\mu}_i)^{\mathsf{T}} $$

c_MLE_mu = np.full_like(data['mu'], np.nan)
c_MLE_cov = np.full_like(data['cov'], np.nan)

for c in range(C):
    c_MLE_mu[c] = np.mean(data['train1'][c], axis=0)
    c_MLE_cov[c] = np.cov(data['train1'][c].T)

show('\mathrm{[Part\ I]\ c)\ \ MLE\ \mu_1^{\mathsf{T}}}', c_MLE_mu[0][:, None].T)
show('\mathrm{[Part\ I]\ c)\ \ MLE\ \mu_2^{\mathsf{T}}}', c_MLE_mu[1][:, None].T)
show('\mathrm{[Part\ I]\ c)\ \ MLE\ \mu_3^{\mathsf{T}}}', c_MLE_mu[2][:, None].T)

$\displaystyle \mathrm{[Part\ I]\ c)\ \ MLE\ \mu_1^{\mathsf{T}}} = \begin{bmatrix} 1.9825 & 2.9415 & 1.3979 & 5.3871 & 9.1868 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ c)\ \ MLE\ \mu_2^{\mathsf{T}}} = \begin{bmatrix} -4.5503 & 5.9388 & -0.6812 & 2.9552 & 10.0478 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ c)\ \ MLE\ \mu_3^{\mathsf{T}}} = \begin{bmatrix} 0.815 & -2.2132 & 1.7862 & -0.4861 & 2.4406 \\ \end{bmatrix}$

show('\mathrm{[Part\ I]\ c)\ \ MLE\ \Sigma_1}', c_MLE_cov[0])
show('\mathrm{[Part\ I]\ c)\ \ MLE\ \Sigma_2}', c_MLE_cov[1])
show('\mathrm{[Part\ I]\ c)\ \ MLE\ \Sigma_3}', c_MLE_cov[2])

$\displaystyle \mathrm{[Part\ I]\ c)\ \ MLE\ \Sigma_1} = \begin{bmatrix} 1.1021 & 0.0069 & 0.2927 & -0.0909 & -0.0012 \\ 0.0069 & 0.5666 & -0.3443 & 0.0503 & 0.0192 \\ 0.2927 & -0.3443 & 2.4534 & -0.1528 & -0.2686 \\ -0.0909 & 0.0503 & -0.1528 & 0.6542 & -0.1083 \\ -1.1985e-03 & 1.9217e-02 & -2.6863e-01 & -1.0834e-01 & 2.3480e+00 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ c)\ \ MLE\ \Sigma_2} = \begin{bmatrix} 2.5306 & -0.4836 & 1.1574 & 0.5788 & 0.0173 \\ -0.4836 & 2.525 & -0.4431 & -0.6565 & 1.058 \\ 1.1574 & -0.4431 & 3.7571 & 1.2025 & 0.3038 \\ 0.5788 & -0.6565 & 1.2025 & 1.5769 & -0.2817 \\ 0.0173 & 1.058 & 0.3038 & -0.2817 & 2.8041 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ c)\ \ MLE\ \Sigma_3} = \begin{bmatrix} 3.2999 & -0.1454 & 0.7844 & 1.8742 & 0.0386 \\ -0.1454 & 7.3366 & 0.2 & -0.1654 & 4.2923 \\ 0.7844 & 0.2 & 3.6297 & 3.3543 & 0.1996 \\ 1.8742 & -0.1654 & 3.3543 & 4.6125 & -0.1246 \\ 0.0386 & 4.2923 & 0.1996 & -0.1246 & 6.443 \\ \end{bmatrix}$

c_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        distances = np.zeros(C)
        for i in range(C):
            distances[i] = discriminant_function(x, c_MLE_mu[i], c_MLE_cov[i])
        pred_class = np.argmax(distances)
        c_confusion_matrix[c][pred_class] += 1

show('\mathrm{[Part\ I]\ c)\ \ Confusion\ matrix}', c_confusion_matrix)

$\displaystyle \mathrm{[Part\ I]\ c)\ \ Confusion\ matrix} = \begin{bmatrix} 9973 & 20 & 7 \\ 6 & 9992 & 2 \\ 21 & 0 & 9979 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ I]\ c)\ accuracy}', accuracy(c_confusion_matrix))
show_percent(r'\mathrm{[Part\ I]\ c)\ error}', 100 - accuracy(c_confusion_matrix))

$\displaystyle \mathrm{[Part\ I]\ c)\ accuracy} = 99.81\%$

$\displaystyle \mathrm{[Part\ I]\ c)\ error} = 0.19\%$

d) Repeat part c), but this time use the Training Data II for the ML and then classify the Testing Data II again.

d_MLE_mu = np.full_like(data['mu'], np.nan)
d_MLE_cov = np.full_like(data['cov'], np.nan)

for c in range(C):
    d_MLE_mu[c] = np.mean(data['train2'][c], axis=0)
    d_MLE_cov[c] = np.cov(data['train2'][c].T)

show('\mathrm{[Part\ I]\ d)\ \ MLE\ \mu_1^{\mathsf{T}}}', d_MLE_mu[0][:, None].T)
show('\mathrm{[Part\ I]\ d)\ \ MLE\ \mu_2^{\mathsf{T}}}', d_MLE_mu[1][:, None].T)
show('\mathrm{[Part\ I]\ d)\ \ MLE\ \mu_3^{\mathsf{T}}}', d_MLE_mu[2][:, None].T)

$\displaystyle \mathrm{[Part\ I]\ d)\ \ MLE\ \mu_1^{\mathsf{T}}} = \begin{bmatrix} 2.0027 & 2.9968 & 1.0015 & 5.4631 & 8.7453 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ d)\ \ MLE\ \mu_2^{\mathsf{T}}} = \begin{bmatrix} -4.417 & 5.9488 & -0.9668 & 2.9721 & 10.0056 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ d)\ \ MLE\ \mu_3^{\mathsf{T}}} = \begin{bmatrix} 1.1914 & -2.2265 & 1.4473 & -0.4403 & 2.5968 \\ \end{bmatrix}$

show('\mathrm{[Part\ I]\ d)\ \ MLE\ \Sigma_1}', d_MLE_cov[0])
show('\mathrm{[Part\ I]\ d)\ \ MLE\ \Sigma_2}', d_MLE_cov[1])
show('\mathrm{[Part\ I]\ d)\ \ MLE\ \Sigma_3}', d_MLE_cov[2])

$\displaystyle \mathrm{[Part\ I]\ d)\ \ MLE\ \Sigma_1} = \begin{bmatrix} 1.0786 & -0.0427 & 0.0524 & 0.0598 & 0.0705 \\ -0.0427 & 0.4941 & -0.0732 & -0.0292 & -0.0506 \\ 0.0524 & -0.0732 & 2.6967 & 0.0661 & 0.1141 \\ 0.0598 & -0.0292 & 0.0661 & 0.6706 & -0.0106 \\ 0.0705 & -0.0506 & 0.1141 & -0.0106 & 3.1491 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ d)\ \ MLE\ \Sigma_2} = \begin{bmatrix} 1.9979 & 0.0135 & 0.8236 & 0.4468 & 0.0927 \\ 0.0135 & 3.5147 & -0.0046 & -0.1368 & 0.3662 \\ 0.8236 & -0.0046 & 4.0323 & 1.0935 & -0.1122 \\ 0.4468 & -0.1368 & 1.0935 & 1.5593 & -0.118 \\ 0.0927 & 0.3662 & -0.1122 & -0.118 & 2.4692 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ I]\ d)\ \ MLE\ \Sigma_3} = \begin{bmatrix} 4.5314 & 0.0381 & 1.6708 & 2.9363 & 1.7751 \\ 0.0381 & 4.711 & -0.1134 & 0.0668 & 3.4286 \\ 1.6708 & -0.1134 & 4.619 & 4.3969 & 0.4896 \\ 2.9363 & 0.0668 & 4.3969 & 5.9432 & 0.5577 \\ 1.7751 & 3.4286 & 0.4896 & 0.5577 & 7.6444 \\ \end{bmatrix}$

d_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        distances = np.zeros(C)
        for i in range(C):
            distances[i] = discriminant_function(x, d_MLE_mu[i], d_MLE_cov[i])
        pred_class = np.argmax(distances)
        d_confusion_matrix[c][pred_class] += 1

show('\mathrm{[Part\ I]\ d)\ \ Confusion\ matrix}', d_confusion_matrix)

$\displaystyle \mathrm{[Part\ I]\ d)\ \ Confusion\ matrix} = \begin{bmatrix} 9981 & 6 & 13 \\ 8 & 9992 & 0 \\ 8 & 0 & 9992 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ I]\ d)\ accuracy}', accuracy(d_confusion_matrix))
show_percent(r'\mathrm{[Part\ I]\ d)\ error}', 100 - accuracy(d_confusion_matrix))

$\displaystyle \mathrm{[Part\ I]\ d)\ accuracy} = 99.88\%$

$\displaystyle \mathrm{[Part\ I]\ d)\ error} = 0.12\%$

e) Make comments on each of the results above and then compare them (’compare’ does not mean to say “this was better than that”, but to say why that was the case).

Part	Question	Algorithm	Training set	Acccuracy on Testing data I	Accuracy on Testing data II
I	a, b	BDR	A prior knowledge	99.87%	99.89%
I	c	BDR	Training data I	-	99.81%
I	d	BDR	Training data II	-	99.88%

Table 2: Part I classification results

First of all, we use a priori knowledge about the data which was given in the instruction, i.e., we know how each sample is drawn from the distribution. In this case, we do not need the training data at all. The accuracy of the questions (a) and (b) are then obviously high (on the two testing sets). When it comes to the questions (c) and (d), instead of using the a priori knowledge, we need to estimate the parameters of the Guassian distributions using maximum likelihood estimation (MLE). Clearly, since the estimated parameters are not the actual parameters, the accuracies on the testing sets should drop. For the question (c) which we train the model on the training data I (smaller number of training samples), the accuracy drops from 99.89% to 99.81%. But if we train the model on the training data II (higher number of training samples), the estimated parameters are more accurate. This results in a higher accuracy of 99.88% (versus the 99.81% accuracy from the training data I). As we all know, the accuracy shall not exceed the accuracy of 99.89% that is from the complete knowledge about the distributions.

---

Part II

In this part, you will reduce the dimensionality of the data by applying MDA. That is,

a) Using MDA on the Training Data I, find the matrix W such that $\vec{y}$ = $W\vec{x}$. (What is the expected dimension of $\vec{y}$?). For each class, compute the $\hat{\vec{\mu}}_i$ and $\Sigma_i$ of the new reduced-space variable $\vec{y}$.

Scatter matrix of class $i$ is defined by: $$ \mathbf{S}_i = \sum_{\vec{x} \in \mathcal{D}_i} (\vec{x} - \vec{m}_i)(\vec{x} - \vec{m}_i)^{\mathsf{T}} $$

The $d$-dimensional sample mean of class $i$ or $\vec{m}_i$ is given by: $$ \vec{m}_i = \frac{1}{n_i} \sum_{\vec{x} \in \mathcal{D}_i} \vec{x} $$ where $n_i$ is the number of training samples in class $i$.

Within-class scatter matrix is the summation of scatter matrices from all classes. $$ \mathbf{S}_W = \sum_{i=1}^{C} \mathbf{S}_i $$

def scatter_matrix(X):
    m = X.mean(axis=0)
    S = np.zeros((len(m), len(m)))
    for x in X:
        S += np.outer(x - m, (x - m).T)
    return S

def within_class_scatter(data):
    S_w = np.zeros((data[0].shape[1], data[0].shape[1]))
    for c in range(len(data)):
        S_w += scatter_matrix(data[c])
    return S_w

Between-class scatter matrix can be computed by: $$ \mathbf{S}_B = \sum_{i=1}^{C} n_i (\vec{m}_i - \vec{m}) (\vec{m}_i - \vec{m})^{\mathsf{T}} $$ where $\vec{m}$ is a total mean vector defined as: $$ \vec{m} = \frac{1}{n} \sum_{\vec{x}} \vec{x} = \frac{1}{n} \sum_{i=1}^{C} n_i \vec{m}_i $$

def total_mean_vector(data):
    m = np.zeros(data[0].shape[1])
    n = 0
    for c in range(len(data)):
        n_i = len(data[c])
        m_i = data[c].mean(axis=0)
        n += n_i
        m += n_i * m_i
    m /= n
    return m

def between_class_scatter(data):
    S_B = np.zeros((data[0].shape[1], data[0].shape[1]))
    m = total_mean_vector(data)
    for c in range(len(data)):
        n_i = len(data[c])
        m_i = data[c].mean(axis=0)
        S_B += n_i * np.outer((m_i - m), (m_i - m).T)
    return S_B

def total_scatter(data):
    return within_class_scatter(data) + between_class_scatter(data)

II_a_S_W = within_class_scatter(data['train1'])

show(r'\mathrm{[Part\ II]\ a)\ \ Within-class\ scatter\ matrix}\ \ S_W', II_a_S_W)

$\displaystyle \mathrm{[Part\ II]\ a)\ \ Within-class\ scatter\ matrix}\ \ S_W = \begin{bmatrix} 339.6975 & -30.4845 & 109.4908 & 115.7368 & 2.6803 \\ -30.4845 & 510.9815 & -28.784 & -37.8109 & 263.1046 \\ 109.4908 & -28.784 & 482.1722 & 215.7956 & 11.5063 \\ 115.7368 & -37.8109 & 215.7956 & 335.3344 & -25.2126 \\ 2.6803 & 263.1046 & 11.5063 & -25.2126 & 568.1629 \\ \end{bmatrix}$

II_a_S_W_inv = np.linalg.inv(II_a_S_W)

show(r'\mathrm{[Part\ II]\ a)\ \ Inverse\ of\ within-class\ scatter\ matrix}\ \ S_W^{-1}', II_a_S_W_inv)

$\displaystyle \mathrm{[Part\ II]\ a)\ \ Inverse\ of\ within-class\ scatter\ matrix}\ \ S_W^{-1} = \begin{bmatrix} 0.0034 & 0.0002 & -0.0003 & -0.0009 & -0.0001 \\ 1.8377e-04 & 2.5991e-03 & 1.1220e-04 & 6.6926e-05 & -1.2037e-03 \\ -0.0003 & 0.0001 & 0.003 & -0.0018 & -0.0002 \\ -9.4591e-04 & 6.6926e-05 & -1.7940e-03 & 4.4863e-03 & 2.0889e-04 \\ -0.0001 & -0.0012 & -0.0002 & 0.0002 & 0.0023 \\ \end{bmatrix}$

II_a_S_B = between_class_scatter(data['train1'])

show(r'\mathrm{[Part\ II]\ a)\ \ Between-class\ scatter\ matrix}\ \ S_B', II_a_S_B)

$\displaystyle \mathrm{[Part\ II]\ a)\ \ Between-class\ scatter\ matrix}\ \ S_B = \begin{bmatrix} 1213.7516 & -955.0083 & 439.46 & 71.3476 & -642.72 \\ -955.0083 & 1700.1563 & -472.4702 & 850.633 & 1656.1391 \\ 439.46 & -472.4702 & 176.0326 & -95.2565 & -386.3353 \\ 71.3476 & 850.633 & -95.2565 & 870.8586 & 1061.7677 \\ -642.72 & 1656.1391 & -386.3353 & 1061.7677 & 1735.3519 \\ \end{bmatrix}$

Solve for the eigenvalues and eigenvectors $$ \mathbf{S_W^{-1}} \mathbf{S_B} \mathbf{w} = \lambda \mathbf{w} $$

II_a_eigenvalues, II_a_eigenvectors = np.linalg.eigh(II_a_S_W_inv.dot(II_a_S_B))

show(r'\mathrm{[Part\ II]\ a)\ \ Eigenvalues}\ \Lambda', II_a_eigenvalues[:, None].T)
show(r'\mathrm{[Part\ II]\ a)\ \ Eigenvectors}\ \Phi', II_a_eigenvectors)

$\displaystyle \mathrm{[Part\ II]\ a)\ \ Eigenvalues}\ \Lambda = \begin{bmatrix} -3.4071 & 0.2572 & 1.6272 & 3.0697 & 11.8115 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ a)\ \ Eigenvectors}\ \Phi = \begin{bmatrix} -0.0274 & 0.0161 & -0.1336 & -0.9447 & -0.2979 \\ 0.7665 & 0.1058 & -0.1324 & -0.1879 & 0.5903 \\ 0.3518 & -0.8832 & 0.1829 & 0.0273 & -0.2488 \\ -0.5232 & -0.4529 & -0.2746 & -0.1582 & 0.6486 \\ -0.1191 & 0.0577 & 0.9251 & -0.2157 & 0.2831 \\ \end{bmatrix}$

For the $C$-class classification problem, we select the largest $C-1$ non-zero eigenvalues to indicate which columns we should use to form a weight matrix $\mathbf{W}$.

II_a_largest_columns = II_a_eigenvalues.argsort()[-(C-1):][::-1]

show(r'\mathrm{[Part\ II]\ a)\ \ Largest}\ C-1 = 3-1 = 2\ \mathrm{eigenvalues, so\ we\ select\ the\ columns}', II_a_largest_columns[:, None].T, precision=0)

$\displaystyle \mathrm{[Part\ II]\ a)\ \ Largest}\ C-1 = 3-1 = 2\ \mathrm{eigenvalues, so\ we\ select\ the\ columns} = \begin{bmatrix} 4 & 3 \\ \end{bmatrix}$

II_a_W = II_a_eigenvectors[:, II_a_largest_columns]
show(r'\mathrm{[Part\ II]\ a)\ \ Weight\ matrix}\ W', II_a_W)

$\displaystyle \mathrm{[Part\ II]\ a)\ \ Weight\ matrix}\ W = \begin{bmatrix} -0.2979 & -0.9447 \\ 0.5903 & -0.1879 \\ -0.2488 & 0.0273 \\ 0.6486 & -0.1582 \\ 0.2831 & -0.2157 \\ \end{bmatrix}$

With the $d$-by-$(C-1)$ weight matrix $\mathbf{W}$, we can project from $d$-dimensional space to a $(C-1)$-dimensional space by taking the dot product between a sample $\vec{x}$ with the weight matrix $\mathbf{W}$. In this case, the expected dimension of $\vec{y}$ after the projection of an $\vec{x}$ with $W$ is $C-1 = 2$ (since $C=3$ is the number of classes for this dataset).

This can also be verified by applying the dot product between a sample $\vec{x}$ and the weight matrix $W$ as follows:

x = data['train1'][0][0]
y = (x @ II_a_W)
print("x's shape is", x[:, None].shape)
print("y's shape is", y[:, None].shape)

x's shape is (5, 1)
y's shape is (2, 1)

The mean vectors and covariance matrices are the estimated just like in the part I, but we estimate on the transformed samples.

II_a_mu = np.full((C, C-1), np.nan)
II_a_cov = np.full((C, C-1, C-1), np.nan)

for c in range(C):
    II_a_Y = data['train1'][c] @ II_a_W
    II_a_mu[c] = II_a_Y.mean(axis=0)
    II_a_cov[c] = np.cov(II_a_Y.T)

show('\mathrm{[Part\ II]\ a)\ \ MLE\ }\mu_1^{\mathsf{T}}', II_a_mu[0][:, None].T)
show('\mathrm{[Part\ II]\ a)\ \ MLE\ }\mu_2^{\mathsf{T}}', II_a_mu[1][:, None].T)
show('\mathrm{[Part\ II]\ a)\ \ MLE\ }\mu_3^{\mathsf{T}}', II_a_mu[2][:, None].T)

$\displaystyle \mathrm{[Part\ II]\ a)\ \ MLE\ }\mu_1^{\mathsf{T}} = \begin{bmatrix} 6.8933 & -5.2215 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ a)\ \ MLE\ }\mu_2^{\mathsf{T}} = \begin{bmatrix} 9.7925 & 0.5289 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ a)\ \ MLE\ }\mu_3^{\mathsf{T}} = \begin{bmatrix} -1.6179 & -0.7549 \\ \end{bmatrix}$

show('\mathrm{[Part\ II]\ a)\ \ MLE\ }\Sigma_1', II_a_cov[0])
show('\mathrm{[Part\ II]\ a)\ \ MLE\ }\Sigma_2', II_a_cov[1])
show('\mathrm{[Part\ II]\ a)\ \ MLE\ }\Sigma_3', II_a_cov[2])

$\displaystyle \mathrm{[Part\ II]\ a)\ \ MLE\ }\Sigma_1 = \begin{bmatrix} 1.1804 & 0.0936 \\ 0.0936 & 1.0958 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ a)\ \ MLE\ }\Sigma_2 = \begin{bmatrix} 1.787 & 0.3096 \\ 0.3096 & 2.4871 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ a)\ \ MLE\ }\Sigma_3 = \begin{bmatrix} 5.0605 & -2.0936 \\ -2.0936 & 4.4021 \\ \end{bmatrix}$

b) Apply the transformation W above to the Testing Data II and classify it using Bayes decision rule. Compute the confusion matrix.

The transformation is simply the matrix-vector product as follows: $$ \vec{y} = \mathbf{W}^{\mathsf{T}} \vec{x} $$

We then classify the transformed sample $\vec{y}$ using the discriminant function

II_a_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        y = x @ II_a_W
        distances = np.zeros(C)
        for i in range(C):
            distances[i] = discriminant_function(y, II_a_mu[i], II_a_cov[i])
        pred_class = np.argmax(distances)
        II_a_confusion_matrix[c][pred_class] += 1

show('\mathrm{[Part\ II]\ b)\ \ Confusion\ matrix}', II_a_confusion_matrix)

$\displaystyle \mathrm{[Part\ II]\ b)\ \ Confusion\ matrix} = \begin{bmatrix} 9943 & 38 & 19 \\ 43 & 9949 & 8 \\ 38 & 0 & 9962 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ II]\ b)\ accuracy}', accuracy(II_a_confusion_matrix))
show_percent(r'\mathrm{[Part\ II]\ b)\ error}', 100 - accuracy(II_a_confusion_matrix))

$\displaystyle \mathrm{[Part\ II]\ b)\ accuracy} = 99.51\%$

$\displaystyle \mathrm{[Part\ II]\ b)\ error} = 0.49\%$

c) Repeat parts a) above, but this time use the Training Data II to compute the matrix W and the $(\hat{\vec{\mu}}_i$ and $\Sigma_i$ for each class of the new reduced-space variable $\vec{y}$. Then classify the Testing Data II again. Compute the confusion matrix.

II_c_S_W = within_class_scatter(data['train2'])
II_c_S_W_inv = np.linalg.inv(II_c_S_W)
II_c_S_B = between_class_scatter(data['train2'])

show(r'\mathrm{[Part\ II]\ c)\ \ Within-class\ scatter\ matrix}\ \ S_W', II_c_S_W)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Within-class\ scatter\ matrix}\ \ S_W = \begin{bmatrix} 3796.3231 & 4.4118 & 1270.878 & 1718.0063 & 967.1779 \\ 4.4118 & 4351.1681 & -95.4054 & -49.5293 & 1868.3428 \\ 1270.878 & -95.4054 & 5662.6528 & 2772.6569 & 245.2722 \\ 1718.0063 & -49.5293 & 2772.6569 & 4078.3687 & 214.1474 \\ 967.1779 & 1868.3428 & 245.2722 & 214.1474 & 6618.0779 \\ \end{bmatrix}$

show(r'\mathrm{[Part\ II]\ c)\ \ Inverse\ of\ within-class\ scatter\ matrix}\ \ S_W^{-1}', II_c_S_W_inv)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Inverse\ of\ within-class\ scatter\ matrix}\ \ S_W^{-1} = \begin{bmatrix} 3.3999e-04 & 1.9670e-05 & -7.6047e-06 & -1.3515e-04 & -5.0584e-05 \\ 1.9670e-05 & 2.6301e-04 & 5.7989e-06 & -4.9817e-06 & -7.7180e-05 \\ -7.6047e-06 & 5.7989e-06 & 2.6514e-04 & -1.7673e-04 & -4.6332e-06 \\ -1.3515e-04 & -4.9817e-06 & -1.7673e-04 & 4.2148e-04 & 1.4070e-05 \\ -5.0584e-05 & -7.7180e-05 & -4.6332e-06 & 1.4070e-05 & 1.8000e-04 \\ \end{bmatrix}$

show(r'\mathrm{[Part\ II]\ c)\ \ Between-class\ scatter\ matrix}\ \ S_B', II_c_S_B)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Between-class\ scatter\ matrix}\ \ S_B = \begin{bmatrix} 12220.5886 & -10093.8603 & 4302.2158 & 273.909 & -7442.2498 \\ -10093.8603 & 17139.0258 & -4645.8919 & 8563.239 & 16067.6857 \\ 4302.2158 & -4645.8919 & 1650.1551 & -994.4278 & -3851.257 \\ 273.909 & 8563.239 & -994.4278 & 8783.3337 & 9739.9227 \\ -7442.2498 & 16067.6857 & -3851.257 & 9739.9227 & 15713.9044 \\ \end{bmatrix}$

II_c_eigenvalues, II_c_eigenvectors = np.linalg.eigh(II_c_S_W_inv.dot(II_c_S_B))
II_c_largest_columns = II_c_eigenvalues.argsort()[-(C-1):][::-1]

show(r'\mathrm{[Part\ II]\ c)\ \ Eigenvalues}\ \Lambda', II_c_eigenvalues[:, None].T)
show(r'\mathrm{[Part\ II]\ c)\ \ Eigenvectors}\ \Phi', II_c_eigenvectors)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Eigenvalues}\ \Lambda = \begin{bmatrix} -3.023 & 0.1398 & 1.5231 & 2.9598 & 12.2893 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Eigenvectors}\ \Phi = \begin{bmatrix} -0.0684 & 0.0146 & 0.1415 & 0.9106 & -0.3819 \\ 0.7381 & -0.0713 & -0.01 & 0.308 & 0.596 \\ 0.3281 & 0.9062 & 0.0193 & -0.0968 & -0.2477 \\ -0.5614 & 0.4027 & -0.2988 & 0.2527 & 0.6078 \\ -0.1664 & 0.106 & 0.9435 & -0.0513 & 0.2612 \\ \end{bmatrix}$

show(r'\mathrm{[Part\ II]\ c)\ \ Largest}\ C-1 = 3-1 = 2\ \mathrm{eigenvalues, so\ we\ select\ the\ columns}', II_c_largest_columns[:, None].T, precision=0)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Largest}\ C-1 = 3-1 = 2\ \mathrm{eigenvalues, so\ we\ select\ the\ columns} = \begin{bmatrix} 4 & 3 \\ \end{bmatrix}$

II_c_W = II_c_eigenvectors[:, II_c_largest_columns]
show(r'\mathrm{[Part\ II]\ c)\ \ Weight\ matrix}\ W', II_c_W)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Weight\ matrix}\ W = \begin{bmatrix} -0.3819 & 0.9106 \\ 0.596 & 0.308 \\ -0.2477 & -0.0968 \\ 0.6078 & 0.2527 \\ 0.2612 & -0.0513 \\ \end{bmatrix}$

With the $d$-by-$(C-1)$ weight matrix $\mathbf{W}$, we can project from $d$-dimensional space to a $(C-1)$-dimensional space by taking the dot product between a sample $\vec{x}$ with the weight matrix $\mathbf{W}$. In this case, the expected dimension of $\vec{y}$ after the projection of an $\vec{x}$ with $W$ is $C-1 = 2$ (since $C=3$ is the number of classes for this dataset).

This can also be verified by applying the dot product between a sample $\vec{x}$ and the weight matrix $W$ as follows:

x = data['train1'][0][0]
y = (x @ II_c_W)
print("x's shape is", x[:, None].shape)
print("y's shape is", y[:, None].shape)

x's shape is (5, 1)
y's shape is (2, 1)

The mean vectors and covariance matrices are the estimated just like in the part I, but we estimate on the transformed samples.

II_c_mu = np.full((C, C-1), np.nan)
II_c_cov = np.full((C, C-1, C-1), np.nan)

for c in range(C):
    Y = data['train2'][c] @ II_c_W
    II_c_mu[c] = Y.mean(axis=0)
    II_c_cov[c] = np.cov(Y.T)

show('\mathrm{[Part\ II]\ c)\ \ MLE\ }\mu_1^{\mathsf{T}}', II_c_mu[0][:, None].T)
show('\mathrm{[Part\ II]\ c)\ \ MLE\ }\mu_2^{\mathsf{T}}', II_c_mu[1][:, None].T)
show('\mathrm{[Part\ II]\ c)\ \ MLE\ }\mu_3^{\mathsf{T}}', II_c_mu[2][:, None].T)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ MLE\ }\mu_1^{\mathsf{T}} = \begin{bmatrix} 6.3778 & 3.5818 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ c)\ \ MLE\ }\mu_2^{\mathsf{T}} = \begin{bmatrix} 9.8912 & -1.8585 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ c)\ \ MLE\ }\mu_3^{\mathsf{T}} = \begin{bmatrix} -1.7297 & 0.0144 \\ \end{bmatrix}$

show('\mathrm{[Part\ II]\ c)\ \ MLE\ }\Sigma_1', II_c_cov[0])
show('\mathrm{[Part\ II]\ c)\ \ MLE\ }\Sigma_2', II_c_cov[1])
show('\mathrm{[Part\ II]\ c)\ \ MLE\ }\Sigma_3', II_c_cov[2])

$\displaystyle \mathrm{[Part\ II]\ c)\ \ MLE\ }\Sigma_1 = \begin{bmatrix} 0.895 & -0.1539 \\ -0.1539 & 1.005 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ c)\ \ MLE\ }\Sigma_2 = \begin{bmatrix} 2.1193 & 0.1692 \\ 0.1692 & 2.1093 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ II]\ c)\ \ MLE\ }\Sigma_3 = \begin{bmatrix} 3.8551 & 1.3983 \\ 1.3983 & 5.2442 \\ \end{bmatrix}$

II_c_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        y = x @ II_c_W
        distances = np.zeros(C)
        for i in range(C):
            distances[i] = discriminant_function(y, II_c_mu[i], II_c_cov[i])
        pred_class = np.argmax(distances)
        II_c_confusion_matrix[c][pred_class] += 1

show('\mathrm{[Part\ II]\ c)\ \ Confusion\ matrix}', II_c_confusion_matrix)

$\displaystyle \mathrm{[Part\ II]\ c)\ \ Confusion\ matrix} = \begin{bmatrix} 9955 & 26 & 19 \\ 48 & 9949 & 3 \\ 26 & 0 & 9974 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ II]\ c)\ \ accuracy}', accuracy(II_c_confusion_matrix))
show_percent(r'\mathrm{[Part\ II]\ c)\ \ error}', 100 - accuracy(II_c_confusion_matrix))

$\displaystyle \mathrm{[Part\ II]\ c)\ \ accuracy} = 99.59\%$

$\displaystyle \mathrm{[Part\ II]\ c)\ \ error} = 0.41\%$

d) Comment on the results from this Part II and then compare these results with the results from the previous Part I.

Part	Question	Algorithm	Training set	Acccuracy on Testing data I	Accuracy on Testing data II
I	a, b	BDR	A prior knowledge	99.87%	99.89%
I	c	BDR	Training data I	-	99.81%
I	d	BDR	Training data II	-	99.88%
II	b	MDA	Training data I	-	99.51%
II	c	MDA	Training data II	-	99.59%

Table 3: Part II classification results

The training data I has smaller number of samples compared to the training data II. Again, this makes the parameter estimation less accurate. The point of Part II is to apply multiple discriminant analysis (MDA), i.e., projection from a $d$-dimensional space to $(C-1)$-dimensional space. This brings about a dimensionality reduction of the feature space. Of course, since the number of features are reduced from 5 to 2 in this case, we would lose some accuracy. However, it is a decent trade-off as the accuracies drop less than 1% when compared to Part I. With the reduced dimensionality, we can apply classifier algorithms more efficient, i.e., it is less computationally expensive.

---

Part III

Now, you will completely forget that you know anything about any of the distributions and/or their parameters and apply a non-parametric approach to classification.

a) First, you will apply Parzen Window to each class in the Testing Data II using the Training Data I and a hypercube window function with $h_n = 0.7$. You must classify the Testing Data II according to the maximum posterior probability. Compute the confusion matrix. Repeat the classification for $h_n = 0.1$ and $h_n = 5$.

Based on the density estimation method, we can estimate an unknown probability density function without knowing its true density and its parameters. A $d$-dimensional hypercube $\mathcal{R}_n$ is a region that we are interested in since we can measure the probability that a sample $\vec{x}$ will fall in the region of class $i$, or $p_n(\vec{n} | \omega_i)$ which is given by: $$ \begin{align*} p_n(\vec{x} | \omega_i) &= \frac{k_n}{n V_n} \end{align*} $$ where $n$ is the total number of samples, and $V_n$ is the volume of the hypercube $\mathcal{R}_n$. So, the volume calculation is quite straightforward according to the definition of the hypercube, we arrive at: $$ V_n = h_n^d $$ where $h_n$ is the length of an edge of the hypercube $\mathcal{R}_n$.

We then define $k_n$ which is the number of samples that reside in a $d$-dimensional hypercube $\mathcal{R}_n$ as: $$ k_n = \sum_{i=1}^{n} \varphi \left( \frac{\vec{x} - \vec{x}_i}{h_n} \right) $$

We combine all of the above equations together, so we could compute the probability by: $$ p_n(\vec{x} | \omega_i) = \frac{k_n}{n V_n} = \frac{1}{n V_n} \sum_{i=1}^{n} \varphi \left( \frac{\vec{x} - \vec{x}_i}{h_n} \right) $$

Based on the decision rule for $C$-class classification problem, a sample $\vec{x}$ is classified to be the predicted class $\hat{\omega}$ by using: $$ \hat{\omega} = \underset{i}{\mathrm{argmax}} p_n(\vec{x} | \omega_i)P(\omega_i) $$ In this project, the priors are equal, thus we can discard the prior term. So, we arrive at the decision rule in terms of class-conditional probability densities: $$ \hat{\omega} = \underset{i}{\mathrm{argmax}} p_n(\vec{x} | \omega_i) $$

$\varphi(\cdot)$ can be any kernel function. The standard hypercube kernel is defined as: $$ \varphi \left(\frac{\vec{x} - \vec{x}_i}{h_n} \right) = \begin{cases} 1 & \text{if}\ |\vec{x} - \vec{x}_i| \leq \frac{h_n}{2} \\ 0 & \text{otherwise} \end{cases} $$

def is_hypercube(kernel_fx):
    return kernel_fx == 'hypercube'

def is_gaussian(kernel_fx):
    return kernel_fx == 'gaussian'

def parzen_window(training_data, x, h_n, kernel_fx='hypercube'):
    d = x.shape[0]
    V_n = h_n ** d
    C = len(training_data)
    if is_hypercube(kernel_fx):
        p_n = np.full(C, np.nan)
    elif is_gaussian(kernel_fx):
        p_n = np.full((C, d), np.nan)
    for c in range(C):
        n = len(training_data[c])
        k = 0
        for x_i in training_data[c]:
            if is_hypercube(kernel_fx) and np.all(np.abs(x - x_i) < h_n / 2):
                k += 1
            elif is_gaussian(kernel_fx):
                k += np.exp(-0.5 * ((x - x_i) / h_n)**2) / ((np.sqrt(2 * np.pi))**d * V_n)
        p_n[c] = (1 / (n * V_n)) * k
    if is_gaussian(kernel_fx):
        p_n = p_n.sum(axis=1)
    pred_class = np.argmax(p_n)
    return pred_class

III_a_01_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train1'], x, h_n=0.1, kernel_fx='hypercube')
        III_a_01_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=0.1: \mathrm{Confusion\ matrix}', III_a_01_confusion_matrix)

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=0.1: \mathrm{Confusion\ matrix} = \begin{bmatrix} 10000 & 0 & 0 \\ 10000 & 0 & 0 \\ 10000 & 0 & 0 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=0.1: \mathrm{accuracy}', accuracy(III_a_01_confusion_matrix))
show_percent(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=0.1: \mathrm{error}', 100 - accuracy(III_a_01_confusion_matrix))

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=0.1: \mathrm{accuracy} = 33.33\%$

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=0.1: \mathrm{error} = 66.67\%$

III_a_07_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train1'], x, h_n=0.7, kernel_fx='hypercube')
        III_a_07_confusion_matrix[c][pred_class] += 1
        
show(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=0.7: \mathrm{Confusion\ matrix}', III_a_07_confusion_matrix)

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=0.7: \mathrm{Confusion\ matrix} = \begin{bmatrix} 10000 & 0 & 0 \\ 9987 & 13 & 0 \\ 9995 & 0 & 5 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=0.7: \mathrm{accuracy}', accuracy(III_a_07_confusion_matrix))
show_percent(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=0.7: \mathrm{error}', 100 - accuracy(III_a_07_confusion_matrix))

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=0.7: \mathrm{accuracy} = 33.39\%$

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=0.7: \mathrm{error} = 66.61\%$

III_a_5_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train1'], x, h_n=5, kernel_fx='hypercube')
        III_a_5_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=5: \mathrm{Confusion\ matrix}', III_a_5_confusion_matrix)

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=5: \mathrm{Confusion\ matrix} = \begin{bmatrix} 9997 & 2 & 1 \\ 442 & 9557 & 1 \\ 1096 & 0 & 8904 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=5: \mathrm{accuracy}', accuracy(III_a_5_confusion_matrix))
show_percent(r'\mathrm{[Part\ III]\ a)\ }\ \ h_n=5: \mathrm{error}', 100 - accuracy(III_a_5_confusion_matrix))

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=5: \mathrm{accuracy} = 94.86\%$

$\displaystyle \mathrm{[Part\ III]\ a)\ }\ \ h_n=5: \mathrm{error} = 5.14\%$

b) Repeat part a), but this time use the Training Data II for the Parzen Window and classify the same Testing Data II with, again, $h_n = 0.1$, $h_n = 0.7$, and $h_n = 5$.

III_b_01_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train2'], x, h_n=0.1, kernel_fx='hypercube')
        III_b_01_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=0.1: \mathrm{Confusion\ matrix}', III_b_01_confusion_matrix)

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=0.1: \mathrm{Confusion\ matrix} = \begin{bmatrix} 10000 & 0 & 0 \\ 10000 & 0 & 0 \\ 10000 & 0 & 0 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=0.1: \mathrm{accuracy}', accuracy(III_b_01_confusion_matrix))
show_percent(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=0.1: \mathrm{error}', 100 - accuracy(III_b_01_confusion_matrix))

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=0.1: \mathrm{accuracy} = 33.33\%$

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=0.1: \mathrm{error} = 66.67\%$

III_b_07_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train2'], x, h_n=0.7, kernel_fx='hypercube')
        III_b_07_confusion_matrix[c][pred_class] += 1
        
show(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=0.7: \mathrm{Confusion\ matrix}', III_b_07_confusion_matrix)

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=0.7: \mathrm{Confusion\ matrix} = \begin{bmatrix} 10000 & 0 & 0 \\ 9819 & 181 & 0 \\ 9924 & 0 & 76 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=0.7: \mathrm{accuracy}', accuracy(III_b_07_confusion_matrix))
show_percent(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=0.7: \mathrm{error}', 100 - accuracy(III_b_07_confusion_matrix))

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=0.7: \mathrm{accuracy} = 34.19\%$

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=0.7: \mathrm{error} = 65.81\%$

III_b_5_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train2'], x, h_n=5, kernel_fx='hypercube')
        III_b_5_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[Part\ III]\ b)\ }\ h_n=5: \mathrm{Confusion\ matrix}', III_b_5_confusion_matrix)

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ h_n=5: \mathrm{Confusion\ matrix} = \begin{bmatrix} 9999 & 0 & 1 \\ 93 & 9906 & 1 \\ 155 & 2 & 9843 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=5: \mathrm{accuracy}', accuracy(III_b_5_confusion_matrix))
show_percent(r'\mathrm{[Part\ III]\ b)\ }\ \ h_n=5: \mathrm{error}', 100 - accuracy(III_b_5_confusion_matrix))

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=5: \mathrm{accuracy} = 99.16\%$

$\displaystyle \mathrm{[Part\ III]\ b)\ }\ \ h_n=5: \mathrm{error} = 0.84\%$

c) Repeat part b) using a Gaussian kernel (a Gaussian window function) with $\sigma = 0.1$. Then repeat this part with $\sigma = 0.7$ and $\sigma = 5$.

In this question, we replace the hypercube kernel with the Gaussian kernel which is given by $$ \varphi \left(\frac{\vec{x} - \vec{x}_i}{\sigma} \right) = \frac{1}{(\sqrt{2\pi})^d V_n} \exp{\left(-\frac{1}{2}(\frac{\vec{x} - \vec{x}_i}{\sigma})^2\right)} $$ where $V_n = \sigma^d$.

III_c_01_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train2'], x, h_n=0.1, kernel_fx='gaussian')
        III_c_01_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.1: \mathrm{Confusion\ matrix}', III_c_01_confusion_matrix)

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.1: \mathrm{Confusion\ matrix} = \begin{bmatrix} 9971 & 27 & 2 \\ 888 & 9109 & 3 \\ 984 & 155 & 8861 \\ \end{bmatrix}$

show_percent(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.1: \mathrm{accuracy}', accuracy(III_c_01_confusion_matrix))
show_percent(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.1: \mathrm{error}', 100 - accuracy(III_c_01_confusion_matrix))

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.1: \mathrm{accuracy} = 93.14\%$

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.1: \mathrm{error} = 6.86\%$

III_c_07_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train2'], x, h_n=0.7, kernel_fx='gaussian')
        III_c_07_confusion_matrix[c][pred_class] += 1
        
show(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.7: \mathrm{Confusion\ matrix}', III_c_07_confusion_matrix)

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.7: \mathrm{Confusion\ matrix} = \begin{bmatrix} 9992 & 8 & 0 \\ 644 & 9355 & 1 \\ 802 & 95 & 9103 \\ \end{bmatrix}$

show_percent(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.7: \mathrm{accuracy}', accuracy(III_c_07_confusion_matrix))
show_percent(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.7: \mathrm{error}', 100 - accuracy(III_c_07_confusion_matrix))

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.7: \mathrm{accuracy} = 94.83\%$

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=0.7: \mathrm{error} = 5.17\%$

III_c_5_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in data['test2'][c]:
        pred_class = parzen_window(data['train2'], x, h_n=5, kernel_fx='gaussian')
        III_c_5_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[III]\ c)\ Gaussian\ kernel} \sigma=5: \mathrm{Confusion\ matrix}', III_c_5_confusion_matrix)

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel} \sigma=5: \mathrm{Confusion\ matrix} = \begin{bmatrix} 10000 & 0 & 0 \\ 219 & 9781 & 0 \\ 745 & 16 & 9239 \\ \end{bmatrix}$

show_percent(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=5: \mathrm{accuracy}', accuracy(III_c_5_confusion_matrix))
show_percent(r'\mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=5: \mathrm{error}', 100 - accuracy(III_c_5_confusion_matrix))

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=5: \mathrm{accuracy} = 96.73\%$

$\displaystyle \mathrm{[III]\ c)\ Gaussian\ kernel}\ \sigma=5: \mathrm{error} = 3.27\%$

d) Comment on the results above. Compare them with the results from the previous Parts above (I and II).

Part	Question	Algorithm	Training set	Acccuracy on Testing data I	Accuracy on Testing data II
I	a, b	BDR	A prior knowledge	99.87%	99.89%
I	c	BDR	Training data I	-	99.81%
I	d	BDR	Training data II	-	99.88%
II	b	MDA	Training data I	-	99.51%
II	c	MDA	Training data II	-	99.59%
III	a	Parzen Window (Hypercube $h_n=0.1$)	Training data I	-	33.33%
III	a	Parzen Window (Hypercube $h_n=0.7$)	Training data I	-	33.39%
III	a	Parzen Window (Hypercube $h_n=5.0$)	Training data I	-	94.86%
III	b	Parzen Window (Hypercube $h_n=0.1$)	Training data II	-	33.33%
III	b	Parzen Window (Hypercube $h_n=0.7$)	Training data II	-	34.19%
III	b	Parzen Window (Hypercube $h_n=5.0$)	Training data II	-	99.16%
III	c	Parzen Window (Gaussian $\sigma=0.1$)	Training data II	-	93.14%
III	c	Parzen Window (Gaussian $\sigma=0.7$)	Training data II	-	94.83%
III	c	Parzen Window (Gaussian $\sigma=5.0$)	Training data II	-	96.73%

Table 4: Part III classification results

The very first setting is that we set the width of the hypercube window to be $h_n = 0.1$ which is way too small. So, there are no samples that are in the Parzen windows at different locations of the training samples. The effect of using argmax in the implementation makes the Parzen window predicts only the class $\omega_1$ which results in only 33.33% accuracy. This is the case for using either training data I or training data II. The results are the same.

When we increase the width of the hypercube kernel to $h_n = 0.7$, the Parzen window luckily covers just some test samples, thus it gains just a little more accuracy (33.39%). The result on the testing data II is 34.19%. The reason is that the number of testing samples are larger. So, we have a little bit more chance to reside in the windows.

But when we set the width of the hypercube kernel to $h_n = 5.0$, the window width is large enough to cover areas that samples might fall in. It yields good classifications results on the test set. It can attain the accuracy of up to 94.86% if we use the training data I. Once we increase the number of training samples (utilizing training data II), it even yields a better accuracy of 99.16%. Due to the fact that we endow the Parzen window with the finite number of training samples, small number of training samples would cause many holes at the fringes (along side with local spikes in the densities). To conclude, the more number of training samples can lead to a better model performance because we have fewer holes at the fringes.

Once we switch over to the Gaussian kernel, even the standard deviation $\sigma=0.1$ produces a very high accuracy of 93.14% due to the characteristic of the Gaussian kernel that has infinite support (unlike the hypercube kernel that has a crisp window). Therefore, even a test sample is far away from the training samples, it still contributes and gives us some value for the term $\varphi \left(\frac{\vec{x} - \vec{x}_i}{h_n} \right)$. So, when it comes to the classification based on the decision rule defined above, we are interested in the class that has the maximum probability, even from infinitesimal values. It seems like if we increase $\sigma$ to 0.7 and 5.0, we get better results, 94.83% and 96.73%, respectively. Because their window widths are larger which make a sample to be more probable to fall in the window of the more likely class.

In terms of the comparison between the Parzen window and the previous classification methods, it seems to be less accurate. However, we need to keep in mind that this is a non-parametric technique which does not require much a priori knowledge. Instead, we need to choose a kernel function and its hyperparameter ($h_n$ or $\sigma$). The biggest problem is that we need to fine-tune or search the best setting for the Parzen window, e.g., we do not know what should be the value of $h_n$. But doing so can be prone to overfitting and generalization issue, we should consider using cross validation to find the best setting. On the flip side, MLE estimates $\vec{\mu}$ and $\Sigma$ of a Gaussian distribution. We see the results from Part I are pretty high since, the assumed distribution matches the true distribution. If the true distribution is not Gaussian and we do not know it, this would be another case. It this the same for Part II which is MDA that provides a dimensionality reduction follwed by a BDR.

Also note that, if we use a large number of training samples, the Parzen window takes more time to execute since its computational complexity involves the number of training samples to go through when performing a classification.

---

Part IV

Once again, you will forget that you know anything about any of the distributions and/or their parameters and apply another non-parametric approach.

a) Using the Training Data I, you will classify the Testing Data II using k-Nearest-Neighbour. Use $k_n = \sqrt{n}$.

Recall that the parzen window method requires us to choose a kernel function (window function) and its size. However, if we do not have that knowledge, we can utilize the k-nearest-neighbor algorithm instead. Unlike Parzen window that fixes the volume $V_n$, k-NN fixes $k$ which is the number of nearest neighbors. This means we grow a cell at a sample $\vec{x}$ until it covers $k_n$ samples. $k_i$ samples out of $k_n$ samples are labeled as class $\omega_i$. So, we can compute the estimated a posteriori probability by: $$ p_n(\omega_i | \vec{x}) = \frac{k_i}{k_n} $$

The way that we find the nearest neighbors is to use a distance function $D$. In this project, euclidean distance which is a Minkowski distance with $p=2$ in $d$-dimensional space is used. $$ D(\vec{a}, \vec{b}) = \left( \sum_{j=1}^{d} (a_j - b_j)^2 \right)^{1/2} $$

Based on the Bayes decision rule, we can classify a sample $\vec{x}$ to be a class $\hat{\omega}$ by computing: $$ \hat{\omega} = \underset{i}{\mathrm{argmax}} p_n(\omega_i | \vec{x}) = \underset{i}{\mathrm{argmax}} k_i $$ which means we classify a sample $\vec{x}$ based on the most frequent class within the grown cell that covers $k_n$ samples.

def squared_euclidean_distance(a, b):
    return np.sum((a - b)**2)

def k_nearest_neighbors(training_data, x, k):
    distances = []
    for c in range(len(training_data)):
        for y in training_data[c]:
            d = squared_euclidean_distance(x, y)
            distances.append((d, c))
    distances = np.array(distances)
    k_indices = np.argpartition(distances[:, 0], k)[:k]
    unique, counts = np.unique(distances[k_indices, 1].astype(int), return_counts=True)
    pred_class = unique[counts.argmax()]
    return pred_class

In this question, since the training data I has 50 samples for each class, this means we have 150 training samples. We can estimate $k_n = \sqrt{n} = \sqrt{150} \approx 12$.

IV_a_training_data = np.array(data['train1'])
IV_a_N = IV_a_training_data.shape[0] * IV_a_training_data.shape[1]
IV_a_k_n = round(np.sqrt(IV_a_N))
show(r'\mathrm{[Part\ IV]\ a)}\ \ k_n', IV_a_k_n, precision=0)

$\displaystyle \mathrm{[Part\ IV]\ a)}\ \ k_n = 12$

IV_a_test_data = np.array(data['test2'])
IV_a_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in IV_a_test_data[c]:
        pred_class = k_nearest_neighbors(IV_a_training_data, x, k=IV_a_k_n)
        IV_a_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[Part\ IV]\ a)}\ \ \mathrm{Confusion\ matrix}', IV_a_confusion_matrix)

$\displaystyle \mathrm{[Part\ IV]\ a)}\ \ \mathrm{Confusion\ matrix} = \begin{bmatrix} 9998 & 0 & 2 \\ 69 & 9930 & 1 \\ 338 & 14 & 9648 \\ \end{bmatrix}$

show_percent(r'\mathrm{[IV]\ a)}\ \ \mathrm{accuracy}', accuracy(IV_a_confusion_matrix))
show_percent(r'\mathrm{[IV]\ a)}\ \ \mathrm{error}', 100 - accuracy(IV_a_confusion_matrix))

$\displaystyle \mathrm{[IV]\ a)}\ \ \mathrm{accuracy} = 98.59\%$

$\displaystyle \mathrm{[IV]\ a)}\ \ \mathrm{error} = 1.41\%$

b) Repeat part a), but this time use the Training Data II to classify the Testing Data II.

In this question, since the training data II has 500 samples for each class, this means we have 1,500 training samples. We can estimate $k_n = \sqrt{n} = \sqrt{1500} \approx 39$.

IV_b_training_data = np.array(data['train2'])
IV_b_N = IV_b_training_data.shape[0] * IV_b_training_data.shape[1]
IV_b_k_n = round(np.sqrt(IV_b_N))
show(r'\mathrm{[Part\ IV]\ b)}\ \ k_n', IV_b_k_n, precision=0)

$\displaystyle \mathrm{[Part\ IV]\ b)}\ \ k_n = 39$

IV_b_test_data = np.array(data['test2'])
IV_b_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in IV_b_test_data[c]:
        pred_class = k_nearest_neighbors(IV_b_training_data, x, k=IV_b_k_n)
        IV_b_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[Part\ IV]\ b)}\ \mathrm{Confusion\ matrix}', IV_b_confusion_matrix)

$\displaystyle \mathrm{[Part\ IV]\ b)}\ \mathrm{Confusion\ matrix} = \begin{bmatrix} 10000 & 0 & 0 \\ 53 & 9946 & 1 \\ 121 & 8 & 9871 \\ \end{bmatrix}$

show_percent(r'\mathrm{[IV]\ b)}\ \ \mathrm{accuracy}', accuracy(IV_b_confusion_matrix))
show_percent(r'\mathrm{[IV]\ b)}\ \ \mathrm{error}', 100 - accuracy(IV_b_confusion_matrix))

$\displaystyle \mathrm{[IV]\ b)}\ \ \mathrm{accuracy} = 99.39\%$

$\displaystyle \mathrm{[IV]\ b)}\ \ \mathrm{error} = 0.61\%$

c) Repeat part b) using your choice for $k_n = f(n)$.

In this question, we define the number of $k_n$ nearest neighbors to be: $$ k_n = f(n) = \frac{\sqrt{n}}{2} $$

Hence, we substitute $n = 1,500$, we arrive at: $$ k_n = \frac{\sqrt{n}}{2} = \frac{\sqrt{1500}}{2} \approx 19 $$

IV_c_k_n = round(np.sqrt(IV_b_N) / 2)
show(r'\mathrm{[Part\ IV]\ c)}\ \ k_n', IV_c_k_n, precision=0)

$\displaystyle \mathrm{[Part\ IV]\ c)}\ \ k_n = 19$

IV_c_confusion_matrix = np.zeros((C, C), dtype='int')

for c in range(C):
    for x in IV_b_test_data[c]:
        pred_class = k_nearest_neighbors(IV_b_training_data, x, k=IV_c_k_n)
        IV_c_confusion_matrix[c][pred_class] += 1

show(r'\mathrm{[Part\ IV]\ c)}\ \mathrm{Confusion\ matrix}', IV_c_confusion_matrix)

$\displaystyle \mathrm{[Part\ IV]\ c)}\ \mathrm{Confusion\ matrix} = \begin{bmatrix} 9998 & 0 & 2 \\ 43 & 9957 & 0 \\ 77 & 5 & 9918 \\ \end{bmatrix}$

show_percent(r'\mathrm{[IV]\ c)}\ \ \mathrm{accuracy}', accuracy(IV_c_confusion_matrix))
show_percent(r'\mathrm{[IV]\ c)}\ \ \mathrm{error}', 100 - accuracy(IV_c_confusion_matrix))

$\displaystyle \mathrm{[IV]\ c)}\ \ \mathrm{accuracy} = 99.58\%$

$\displaystyle \mathrm{[IV]\ c)}\ \ \mathrm{error} = 0.42\%$

d) Comment on the results above. Compare them with the results from the previous Parts above, in special Part III.

Part	Question	Algorithm	Training set	Acccuracy on Testing data I	Accuracy on Testing data II
I	a, b	BDR	A prior knowledge	99.87%	99.89%
I	c	BDR	Training data I	-	99.81%
I	d	BDR	Training data II	-	99.88%
II	b	MDA	Training data I	-	99.51%
II	c	MDA	Training data II	-	99.59%
III	a	Parzen Window (Hypercube $h_n=0.1$)	Training data I	-	33.33%
III	a	Parzen Window (Hypercube $h_n=0.7$)	Training data I	-	33.39%
III	a	Parzen Window (Hypercube $h_n=5.0$)	Training data I	-	94.86%
III	b	Parzen Window (Hypercube $h_n=0.1$)	Training data II	-	33.33%
III	b	Parzen Window (Hypercube $h_n=0.7$)	Training data II	-	34.19%
III	b	Parzen Window (Hypercube $h_n=5.0$)	Training data II	-	99.16%
III	c	Parzen Window (Gaussian $\sigma=0.1$)	Training data II	-	93.14%
III	c	Parzen Window (Gaussian $\sigma=0.7$)	Training data II	-	94.83%
III	c	Parzen Window (Gaussian $\sigma=5.0$)	Training data II	-	96.73%
IV	a	$k_n$-NN ($k_n = \sqrt{n}$)	Training data I	-	98.59%
IV	b	$k_n$-NN ($k_n = \sqrt{n}$)	Training data II	-	99.39%
IV	c	$k_n$-NN ($k_n = \frac{\sqrt{n}}{2}$)	Training data II	-	99.58%

Table 5: Part IV classification results

First of all, we compute $k_n = \sqrt{n} = \sqrt{150} = 12$. As can be seen from the above results, $12$-NN yields 98.59% accuracy when we use the training data I. The accuracy raises to 99.39% when we increase the number of training samples, i.e., using training data II. Therefore, we have 1,500 samples which means we use $39$-NN in this case. Once we change $k_n$ to 19 based on our owned defined function ($k_n = \frac{\sqrt(n)}{2}$), the accuracy increases to 99.58%. All of these literally depend on the choice of $k_n$ that needs empirical experiments to find the best $k_n$. But we need to be careful about the overfitting and generalization problems as a low training error does not guarantee a small test error. We can also do cross validation to mitigate this issue.

Interestingly, if we compare the results from $k_n$-NN and the Parzen window, all the results from $k_n$-NN are better than those from the Parzen windows. This can lead to the fact that, when it comes to the window function and its hyperparameter of the Parzen window, they are unknown in this case. Hence, $k_n$-NN is a better solution compared to the Parzen window. Because $k_n$-nearest neighbors fixes the number of nearest neighbors $k_n$, as opposed to the Parzen window that fixes the volume $V$ (i.e., it requires the window and the width). To conclude, it seems to provide a better classification results without knowing the underlying knowledge about the data (e.g., the underlying distributions). The performance of an $k_n$-NN will be increase as the number of training samples goes up. However, also note that, if we use a large number of training samples, it takes more time to run since its computational complexity involves the number of training samples to go through when performing a classification.

---

Part V

Once again, you will forget that you know anything about any of the distributions and/or their parameters and apply another non-parametric approach. This time, repeat Part IV above using a linear classifier and the Perceptron criterion – for part a) and c), you obviously cannot pick a value for $k_n$, so, instead, use $\eta = \frac{1}{2}$ for part a) and then use $\eta = \frac{1}{\sqrt{k}}$ for part c), where $k$ is the iteration.

a) Using the Training Data I, you will classify the Testing Data II using the Percetron criterion. Use $\eta = \frac{1}{2}$.

V_a_training_data = data['train1']
V_a_test_data = data['test2']
d = len(V_a_training_data[0][0])

We classify a sample $\vec{x}$ based on the linear discriminant function $g_i(\vec{x})$ of class $\omega_i$ which is given by $$ g_i(\vec{x}) = \vec{a}_i^{\mathsf{T}}\vec{y} $$ where $\vec{a}$ is a weight vector, and $\vec{y}$ is a sample or feature vector.

Therefore, our decision rule to assign $\vec{x}$ to the class $\omega_i$ if $g_i(\vec{x}) > g_j(\vec{x})$, $\forall j \neq i$.

The Perceptron Criterion function is defined as: $$ J(\vec{a}) = \sum_{\vec{y} \in \mathcal{Y}} \left( - \vec{a}^{\mathsf{T}} \vec{y} \right) $$ where $\mathcal{Y}$ is the set of misclassified samples.

We can derive the gradient of $J$ w.r.t. the weight vector $\vec{a}$ at iteration $k$ by: $$ \nabla J(\vec{a}(k)) = \frac{\partial J(\vec{a}(k))}{\partial \vec{a}(k)} = \sum_{\vec{y} \in \mathcal{Y}} (- \vec{y}) $$

def perceptron_criterion(Y):
    Y = np.array(Y)
    return (-Y).sum(axis=0)

Since the question does not specify the way to initialize weight vectors, we randomly initialize $C$ weight vectors at time step 0 by drawing from the Gaussian distribution $\vec{a}(0) \sim \mathcal{N}(\vec{\mu},\,\mathbf{\Sigma})\,$ where $\vec{\mu} = \vec{0}$ and $\mathbf{\Sigma} = \mathbf{I}$.

Note that the random seed is fixed, so that the result is reproducible.

np.random.seed(7720)
V_a_weights = np.random.normal(0, 1, (C, d))

show(r'\mathrm{[Part\ V]\ a)\ \ initial\ weight\ vector\ of\ class 1\ :\ }\vec{a}_1^{\mathsf{T}}(0)', V_a_weights[0][:, None].T)
show(r'\mathrm{[Part\ V]\ a)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_2^{\mathsf{T}}(0)', V_a_weights[1][:, None].T)
show(r'\mathrm{[Part\ V]\ a)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_3^{\mathsf{T}}(0)', V_a_weights[2][:, None].T)

$\displaystyle \mathrm{[Part\ V]\ a)\ \ initial\ weight\ vector\ of\ class 1\ :\ }\vec{a}_1^{\mathsf{T}}(0) = \begin{bmatrix} 0.9566 & 0.7696 & 0.9475 & -0.7427 & -0.5379 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ a)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_2^{\mathsf{T}}(0) = \begin{bmatrix} 1.0558 & 1.4355 & 0.2965 & -0.2638 & -0.2981 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ a)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_3^{\mathsf{T}}(0) = \begin{bmatrix} -0.5351 & 0.1132 & 0.1097 & -0.9897 & 2.8318 \\ \end{bmatrix}$

In terms of training the weight vectors, the update rule is employed based on the delta rule and the above gradient. $$ \begin{align*} \vec{a}(k+1) &= \vec{a}(k) - \eta(k) \nabla J(\vec{a}(k))\\ &= \vec{a}(k) - \eta(k) \sum_{\vec{y} \in \mathcal{Y}} (- \vec{y})\\ &= \vec{a}(k) + \eta(k) \sum_{\vec{y} \in \mathcal{Y}} \vec{y} \end{align*} $$

In this miniproject, the stopping criterion for training the Perceptron criterion function is the maximum iteration $K=10$. This means we train the weights of the model 10 iterations and then stop.

def update_weights(a, learning_rate, gradient):
    a -= learning_rate * gradient
    return a

The learning rate $\eta$ is given by the question (a) to be $\eta = \frac{1}{2}$.

V_a_eta = 1/2
V_a_max_epoch = 25

def perceptron_train(weights, training_data, eta, max_epoch, bias=False):
    C = len(weights)
    for epoch in range(1, max_epoch+1):
#     print(f'[Part V] training epoch {epoch}')
        for c in range(C):
            Y = []
            for x in training_data[c]:
                if bias:
                    x = np.concatenate(([1], x))
                linear_discriminant_functions = np.full(C, np.nan)
                for i in range(C):
                    linear_discriminant_functions[i] = x @ weights[i]
                pred_class = linear_discriminant_functions.argmax()
                if pred_class != c:
                    Y.append(x)
            gradient = perceptron_criterion(Y)
            eta = 1 / np.sqrt(epoch) if eta is None else eta
            weights[c] = update_weights(weights[c], eta, gradient)
    return weights

V_a_weights = perceptron_train(V_a_weights, V_a_training_data, V_a_eta, max_epoch=V_a_max_epoch)

show(r'\mathrm{[Part\ V]\ a)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 1}\ : \vec{a}_1(' + str(V_a_max_epoch) +')', V_a_weights[0][:, None].T)
show(r'\mathrm{[Part\ V]\ a)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 2}\ : \vec{a}_2(' + str(V_a_max_epoch) +')', V_a_weights[1][:, None].T)
show(r'\mathrm{[Part\ V]\ a)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 3}\ : \vec{a}_3(' + str(V_a_max_epoch) +')', V_a_weights[2][:, None].T)

$\displaystyle \mathrm{[Part\ V]\ a)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 1}\ : \vec{a}_1(25) = \begin{bmatrix} 54.5929 & 86.7178 & 39.5987 & 153.9011 & 269.0853 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ a)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 2}\ : \vec{a}_2(25) = \begin{bmatrix} -113.7491 & 151.7652 & -15.5262 & 75.5623 & 257.0636 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ a)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 3}\ : \vec{a}_3(25) = \begin{bmatrix} 66.4962 & -4.2774 & 143.7904 & -7.5259 & 307.85 \\ \end{bmatrix}$

def perceptron_predict(weights, data, bias=False):
    C = len(weights)
    cm = np.zeros((C, C), dtype='int')
    for c in range(C):
        Y = []
        for x in data[c]:
            if bias:
                x = np.concatenate(([1], x))
            linear_discriminant_functions = np.full(C, np.nan)
            for i in range(C):
                linear_discriminant_functions[i] = x @ weights[i]
            pred_class = linear_discriminant_functions.argmax()
            if pred_class != c:
                Y.append(x)
            cm[c][pred_class] += 1
    return cm

V_a_training_confusion_matrix = perceptron_predict(V_a_weights, V_a_training_data)

show(r'\mathrm{[Part\ V]\ a)}\ \mathrm{Training\ Confusion\ matrix}', V_a_training_confusion_matrix)

$\displaystyle \mathrm{[Part\ V]\ a)}\ \mathrm{Training\ Confusion\ matrix} = \begin{bmatrix} 50 & 0 & 0 \\ 0 & 50 & 0 \\ 0 & 0 & 50 \\ \end{bmatrix}$

V_a_confusion_matrix = perceptron_predict(V_a_weights, V_a_test_data)

show(r'\mathrm{[Part\ V]\ a)}\ \mathrm{Test\ Confusion\ matrix}', V_a_confusion_matrix)

$\displaystyle \mathrm{[Part\ V]\ a)}\ \mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9988 & 0 & 12 \\ 113 & 9882 & 5 \\ 185 & 290 & 9525 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ V]\ a)}\ \ \mathrm{Test\ accuracy}', accuracy(V_a_confusion_matrix))
show_percent(r'\mathrm{[Part\ V]\ a)}\ \ \mathrm{Test\ error}', 100 - accuracy(V_a_confusion_matrix))

$\displaystyle \mathrm{[Part\ V]\ a)}\ \ \mathrm{Test\ accuracy} = 97.98\%$

$\displaystyle \mathrm{[Part\ V]\ a)}\ \ \mathrm{Test\ error} = 2.02\%$

Let us try including the bias $w_0$ to the weight vector $\vec{a}$.

According to the defined linear discriminant function above, $g_i(\vec{x}) = \vec{a}_i^{\mathsf{T}}\vec{y}$, in order to include a bias term $w_0$ to the weights, our weight vector $\vec{a}$ is modified to be: $$ \vec{a} = \begin{bmatrix} w_0\\ w_1\\ w_2\\ w_3\\ w_4\\ w_5 \end{bmatrix} $$ , and the feature vector becomes an augmented feature vector which is given by:

$$ \vec{y} = \begin{bmatrix} 1\\ x_1\\ x_2\\ x_3\\ x_4\\ x_5 \end{bmatrix} = \begin{bmatrix} w_0\\ \\ \vec{x}\\ \\ \end{bmatrix} $$

Everything else is the same.

np.random.seed(7720)
V_a_weights_with_bias = np.random.normal(0, 1, (C, d+1))
V_a_weights_with_bias = perceptron_train(V_a_weights_with_bias, V_a_training_data, V_a_eta, max_epoch=V_a_max_epoch, bias=True)
V_a_confusion_matrix_with_bias = perceptron_predict(V_a_weights_with_bias, V_a_test_data, bias=True)

show(r'\mathrm{[Part\ V]\ a)\ \ with\ bias}\ \ w_0\ ,K=' + str(V_a_max_epoch) + ':\ \mathrm{Test\ Confusion\ matrix}', V_a_confusion_matrix_with_bias)
show_percent(r'\mathrm{[Part\ V]\ a)\ \ with\ bias}\ \ w_0\ ,K=' + str(V_a_max_epoch) + ':\mathrm{Test\ accuracy}', accuracy(V_a_confusion_matrix_with_bias))

$\displaystyle \mathrm{[Part\ V]\ a)\ \ with\ bias}\ \ w_0\ ,K=25:\ \mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9989 & 0 & 11 \\ 102 & 9895 & 3 \\ 148 & 244 & 9608 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ a)\ \ with\ bias}\ \ w_0\ ,K=25:\mathrm{Test\ accuracy} = 98.31\%$

b) Repeat part a), but this time use the Training Data II to classify the Testing Data II.

np.random.seed(7720)
V_b_weights = np.random.normal(0, 1, (C, d))
V_b_max_epoch = 25

show(r'\mathrm{[Part\ V]\ b)\ \ initial\ weight\ vector\ of\ class 1\ :\ }\vec{a}_1^{\mathsf{T}}(0)', V_b_weights[0][:, None].T)
show(r'\mathrm{[Part\ V]\ b)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_2^{\mathsf{T}}(0)', V_b_weights[1][:, None].T)
show(r'\mathrm{[Part\ V]\ b)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_3^{\mathsf{T}}(0)', V_b_weights[2][:, None].T)

$\displaystyle \mathrm{[Part\ V]\ b)\ \ initial\ weight\ vector\ of\ class 1\ :\ }\vec{a}_1^{\mathsf{T}}(0) = \begin{bmatrix} 0.9566 & 0.7696 & 0.9475 & -0.7427 & -0.5379 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ b)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_2^{\mathsf{T}}(0) = \begin{bmatrix} 1.0558 & 1.4355 & 0.2965 & -0.2638 & -0.2981 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ b)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_3^{\mathsf{T}}(0) = \begin{bmatrix} -0.5351 & 0.1132 & 0.1097 & -0.9897 & 2.8318 \\ \end{bmatrix}$

The learning rate $\eta$ is given by the question (b) to be $\eta = \frac{1}{2}$ (same as the question (a)).

V_b_training_data = data['train2']
V_b_test_data = data['test2']
d = len(V_b_training_data[0][0])
V_b_eta = 1/2

V_b_weights = perceptron_train(V_b_weights, V_b_training_data, V_b_eta, max_epoch=V_b_max_epoch)

V_b_training_confusion_matrix = perceptron_predict(V_b_weights, V_b_training_data)

show(r'\mathrm{[Part\ V]\ b)}\ \mathrm{Training\ Confusion\ matrix}', V_b_training_confusion_matrix)

$\displaystyle \mathrm{[Part\ V]\ b)}\ \mathrm{Training\ Confusion\ matrix} = \begin{bmatrix} 496 & 0 & 4 \\ 2 & 498 & 0 \\ 11 & 18 & 471 \\ \end{bmatrix}$

show(r'\mathrm{[Part\ V]\ b)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 1}\ : \vec{a}_1(' + str(V_b_max_epoch) +')', V_b_weights[0][:, None].T)
show(r'\mathrm{[Part\ V]\ b)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 2}\ : \vec{a}_2(' + str(V_b_max_epoch) +')', V_b_weights[1][:, None].T)
show(r'\mathrm{[Part\ V]\ b)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 3}\ : \vec{a}_3(' + str(V_b_max_epoch) +')', V_b_weights[2][:, None].T)

$\displaystyle \mathrm{[Part\ V]\ b)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 1}\ : \vec{a}_1(25) = \begin{bmatrix} 606.2316 & 951.4675 & 379.4202 & 1700.1509 & 2892.9734 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ b)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 2}\ : \vec{a}_2(25) = \begin{bmatrix} -1153.4315 & 1593.3298 & -163.7283 & 880.7507 & 2811.5055 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ b)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 3}\ : \vec{a}_3(25) = \begin{bmatrix} 898.8788 & -314.1971 & 1364.898 & -147.7004 & 3529.7146 \\ \end{bmatrix}$

V_b_confusion_matrix = perceptron_predict(V_b_weights, V_b_test_data)

show(r'\mathrm{[Part\ V]\ b)}\ \mathrm{Test\ Confusion\ matrix}', V_b_confusion_matrix)

$\displaystyle \mathrm{[Part\ V]\ b)}\ \mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9951 & 0 & 49 \\ 91 & 9904 & 5 \\ 154 & 229 & 9617 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ V]\ b)}\ \ \mathrm{Test\ accuracy}', accuracy(V_b_confusion_matrix))
show_percent(r'\mathrm{[Part\ V]\ b)}\ \ \mathrm{Test\ error}', 100 - accuracy(V_b_confusion_matrix))

$\displaystyle \mathrm{[Part\ V]\ b)}\ \ \mathrm{Test\ accuracy} = 98.24\%$

$\displaystyle \mathrm{[Part\ V]\ b)}\ \ \mathrm{Test\ error} = 1.76\%$

np.random.seed(7720)
V_b_weights_with_bias = np.random.normal(0, 1, (C, d+1))
V_b_weights_with_bias = perceptron_train(V_b_weights_with_bias, V_b_training_data, V_b_eta, max_epoch=25, bias=True)
V_b_confusion_matrix_with_bias = perceptron_predict(V_b_weights_with_bias, V_b_test_data, bias=True)

show(r'\mathrm{[Part\ V]\ b)\ \ with\ bias}\ \ w_0\ ,K=' + str(V_b_max_epoch) + ':\ \mathrm{Test\ Confusion\ matrix}', V_b_confusion_matrix_with_bias)
show_percent(r'\mathrm{[Part\ V]\ b)\ \ with\ bias}\ \ w_0\ ,K=' + str(V_b_max_epoch) + ':\mathrm{Test\ accuracy}', accuracy(V_b_confusion_matrix_with_bias))

$\displaystyle \mathrm{[Part\ V]\ b)\ \ with\ bias}\ \ w_0\ ,K=25:\ \mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9992 & 1 & 7 \\ 75 & 9925 & 0 \\ 117 & 292 & 9591 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ b)\ \ with\ bias}\ \ w_0\ ,K=25:\mathrm{Test\ accuracy} = 98.36\%$

c) Repeat part b) using $\eta = \frac{1}{\sqrt{k}}$.

np.random.seed(7720)
V_c_weights = np.random.uniform(0, 1, (C, d))
V_c_max_epoch = 25

show(r'\mathrm{[Part\ V]\ c)\ \ initial\ weight\ vector\ of\ class 1\ :\ }\vec{a}_1^{\mathsf{T}}(0)', V_c_weights[0][:, None].T)
show(r'\mathrm{[Part\ V]\ c)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_2^{\mathsf{T}}(0)', V_c_weights[1][:, None].T)
show(r'\mathrm{[Part\ V]\ c)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_3^{\mathsf{T}}(0)', V_c_weights[2][:, None].T)

$\displaystyle \mathrm{[Part\ V]\ c)\ \ initial\ weight\ vector\ of\ class 1\ :\ }\vec{a}_1^{\mathsf{T}}(0) = \begin{bmatrix} 0.715 & 0.7672 & 0.2853 & 0.7739 & 0.8136 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_2^{\mathsf{T}}(0) = \begin{bmatrix} 0.3402 & 0.5591 & 0.7862 & 0.1401 & 0.1815 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ initial\ weight\ vector\ of\ class 2\ :\ }\vec{a}_3^{\mathsf{T}}(0) = \begin{bmatrix} 0.596 & 0.0461 & 0.0861 & 0.0813 & 0.1122 \\ \end{bmatrix}$

V_c_training_data = data['train2']
V_c_test_data = data['test2']
d = len(V_c_training_data[0][0])

The learning rate $\eta$ is given by the question (c) to be $\eta = \frac{1}{\sqrt{k}}$.

V_c_weights = perceptron_train(V_c_weights, V_c_training_data, eta=None, max_epoch=V_c_max_epoch)

show(r'\mathrm{[Part\ V]\ c)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 1}\ : \vec{a}_1(' + str(V_c_max_epoch) +')', V_c_weights[0][:, None].T)
show(r'\mathrm{[Part\ V]\ c)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 2}\ : \vec{a}_2(' + str(V_c_max_epoch) +')', V_c_weights[1][:, None].T)
show(r'\mathrm{[Part\ V]\ c)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 3}\ : \vec{a}_3(' + str(V_c_max_epoch) +')', V_c_weights[2][:, None].T)

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 1}\ : \vec{a}_1(25) = \begin{bmatrix} 1214.92 & 1901.833 & 761.1689 & 3402.4621 & 5791.5207 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 2}\ : \vec{a}_2(25) = \begin{bmatrix} -2310.9433 & 3191.131 & -326.6182 & 1767.9339 & 5633.0903 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ \ \mathrm{trained\ weight\ vector\ of\ class\ 3}\ : \vec{a}_3(25) = \begin{bmatrix} 1760.5867 & -530.8207 & 2730.62 & -262.9001 & 7048.0308 \\ \end{bmatrix}$

V_c_training_confusion_matrix = perceptron_predict(V_c_weights, V_c_training_data)

show(r'\mathrm{[Part\ V]\ c)}\ \mathrm{Training\ Confusion\ matrix}', V_c_training_confusion_matrix)

$\displaystyle \mathrm{[Part\ V]\ c)}\ \mathrm{Training\ Confusion\ matrix} = \begin{bmatrix} 496 & 0 & 4 \\ 2 & 498 & 0 \\ 11 & 19 & 470 \\ \end{bmatrix}$

V_c_confusion_matrix = perceptron_predict(V_c_weights, V_c_test_data)

show(r'\mathrm{[Part\ V]\ c)}\ \mathrm{Test\ Confusion\ matrix}', V_c_confusion_matrix)

$\displaystyle \mathrm{[Part\ V]\ c)}\ \mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9946 & 0 & 54 \\ 87 & 9908 & 5 \\ 158 & 231 & 9611 \\ \end{bmatrix}$

show_percent(r'\mathrm{[Part\ V]\ c)}\ \ K=' + str(V_c_max_epoch) +'\ :\mathrm{Test\ accuracy}', accuracy(V_c_confusion_matrix))
show_percent(r'\mathrm{[Part\ V]\ c)}\ \ K=' + str(V_c_max_epoch) +'\ :\mathrm{Test\ error}', 100 - accuracy(V_c_confusion_matrix))

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ K=25\ :\mathrm{Test\ accuracy} = 98.22\%$

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ K=25\ :\mathrm{Test\ error} = 1.78\%$

np.random.seed(7720)
V_c_weights_with_bias = np.random.normal(0, 1, (C, d+1))
V_c_weights_with_bias = perceptron_train(V_c_weights_with_bias, V_c_training_data, eta=None, max_epoch=V_c_max_epoch, bias=True)
V_c_confusion_matrix_with_bias = perceptron_predict(V_c_weights_with_bias, V_c_test_data, bias=True)

show(r'\mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=' + str(V_c_max_epoch) + ':\mathrm{Test\ Confusion\ matrix}', V_c_confusion_matrix_with_bias)
show_percent(r'\mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=' + str(V_c_max_epoch) + ':\mathrm{Test\ accuracy}', accuracy(V_c_confusion_matrix_with_bias))

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=25:\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9993 & 0 & 7 \\ 75 & 9925 & 0 \\ 116 & 292 & 9592 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=25:\mathrm{Test\ accuracy} = 98.37\%$

Let us try to change our only hyperparameter here, the maximum number of iterations $K \in \{10, 20, 40, 60, 80\}$.

for max_iteration in (10, 20, 40, 60, 80):
    np.random.seed(7720)
    V_c_vary_weights = np.random.uniform(0, 1, (C, d))
    V_c_vary_weights = perceptron_train(V_c_vary_weights, V_c_training_data, eta=None, max_epoch=max_iteration)
    V_c_vary_confusion_matrix = perceptron_predict(V_c_vary_weights, V_c_test_data)

    show(r'\mathrm{[Part\ V]\ c)}\ \ K=' + str(max_iteration) + '\ :\mathrm{Test\ Confusion\ matrix}', V_c_vary_confusion_matrix)
    show_percent(r'\mathrm{[V]\ c)}\ \ K=' + str(max_iteration) + '\ :\mathrm{Test\ accuracy}', accuracy(V_c_vary_confusion_matrix))

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ K=10\ :\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9995 & 2 & 3 \\ 64 & 9936 & 0 \\ 140 & 389 & 9471 \\ \end{bmatrix}$

$\displaystyle \mathrm{[V]\ c)}\ \ K=10\ :\mathrm{Test\ accuracy} = 98.01\%$

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ K=20\ :\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9954 & 0 & 46 \\ 85 & 9914 & 1 \\ 141 & 257 & 9602 \\ \end{bmatrix}$

$\displaystyle \mathrm{[V]\ c)}\ \ K=20\ :\mathrm{Test\ accuracy} = 98.23\%$

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ K=40\ :\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9940 & 3 & 57 \\ 92 & 9901 & 7 \\ 180 & 196 & 9624 \\ \end{bmatrix}$

$\displaystyle \mathrm{[V]\ c)}\ \ K=40\ :\mathrm{Test\ accuracy} = 98.22\%$

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ K=60\ :\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9957 & 4 & 39 \\ 89 & 9904 & 7 \\ 198 & 188 & 9614 \\ \end{bmatrix}$

$\displaystyle \mathrm{[V]\ c)}\ \ K=60\ :\mathrm{Test\ accuracy} = 98.25\%$

$\displaystyle \mathrm{[Part\ V]\ c)}\ \ K=80\ :\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9957 & 6 & 37 \\ 92 & 9899 & 9 \\ 207 & 182 & 9611 \\ \end{bmatrix}$

$\displaystyle \mathrm{[V]\ c)}\ \ K=80\ :\mathrm{Test\ accuracy} = 98.22\%$

Let us try varing the maximum number of iterations $K \in \{10, 20, 40, 60, 80\}$ for the models with the bias term.

for max_iteration in (10, 20, 40, 60, 80):
    np.random.seed(7720)
    V_c_vary_weights_with_bias = np.random.uniform(0, 1, (C, d+1))
    V_c_vary_weights_with_bias = perceptron_train(V_c_vary_weights_with_bias, V_c_training_data, eta=None, max_epoch=max_iteration, bias=True)
    V_c_vary_confusion_matrix_with_bias = perceptron_predict(V_c_vary_weights_with_bias, V_c_test_data, bias=True)

    show(r'\mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=' + str(max_iteration) + ':\mathrm{Test\ Confusion\ matrix}', V_c_vary_confusion_matrix_with_bias)
    show_percent(r'\mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=' + str(max_iteration) + ':\mathrm{Test\ accuracy}', accuracy(V_c_vary_confusion_matrix_with_bias))

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=10:\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9986 & 11 & 3 \\ 41 & 9959 & 0 \\ 128 & 348 & 9524 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=10:\mathrm{Test\ accuracy} = 98.23\%$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=20:\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9963 & 0 & 37 \\ 86 & 9912 & 2 \\ 119 & 218 & 9663 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=20:\mathrm{Test\ accuracy} = 98.46\%$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=40:\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9960 & 2 & 38 \\ 91 & 9903 & 6 \\ 134 & 166 & 9700 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=40:\mathrm{Test\ accuracy} = 98.54\%$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=60:\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9966 & 4 & 30 \\ 90 & 9904 & 6 \\ 139 & 152 & 9709 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=60:\mathrm{Test\ accuracy} = 98.60\%$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=80:\mathrm{Test\ Confusion\ matrix} = \begin{bmatrix} 9970 & 4 & 26 \\ 89 & 9904 & 7 \\ 150 & 140 & 9710 \\ \end{bmatrix}$

$\displaystyle \mathrm{[Part\ V]\ c)\ \ with\ bias}\ \ w_0\ ,K=80:\mathrm{Test\ accuracy} = 98.61\%$

d) Comment on the results above. Compare them with the results from the previous Parts above.

Part	Question	Algorithm	Training set	Acccuracy on Testing data I	Accuracy on Testing data II
I	a, b	BDR	A prior knowledge	99.87%	99.89%
I	c	BDR	Training data I	-	99.81%
I	d	BDR	Training data II	-	99.88%
II	b	MDA	Training data I	-	99.51%
II	c	MDA	Training data II	-	99.59%
III	a	Parzen Window (Hypercube $h_n=0.1$)	Training data I	-	33.33%
III	a	Parzen Window (Hypercube $h_n=0.7$)	Training data I	-	33.39%
III	a	Parzen Window (Hypercube $h_n=5.0$)	Training data I	-	94.86%
III	b	Parzen Window (Hypercube $h_n=0.1$)	Training data II	-	33.33%
III	b	Parzen Window (Hypercube $h_n=0.7$)	Training data II	-	34.19%
III	b	Parzen Window (Hypercube $h_n=5.0$)	Training data II	-	99.16%
III	c	Parzen Window (Gaussian $\sigma=0.1$)	Training data II	-	93.14%
III	c	Parzen Window (Gaussian $\sigma=0.7$)	Training data II	-	94.83%
III	c	Parzen Window (Gaussian $\sigma=5.0$)	Training data II	-	96.73%
IV	a	$k_n$-NN ($k_n = \sqrt{n}$)	Training data I	-	98.59%
IV	b	$k_n$-NN ($k_n = \sqrt{n}$)	Training data II	-	99.39%
IV	c	$k_n$-NN ($k_n = \frac{\sqrt{n}}{2}$)	Training data II	-	99.58%
V	a	The Perceptron Criterion ($\eta=\frac{1}{2}$), $K=25$	Training data I	-	97.98%
V	b	The Perceptron Criterion ($\eta=\frac{1}{2}$), $K=25$	Training data II	-	98.24%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=25$	Training data II	-	98.22%
V	a	The Perceptron Criterion ($\eta=\frac{1}{2}$), $K=25$, with bias	Training data I	-	98.31%
V	b	The Perceptron Criterion ($\eta=\frac{1}{2}$), $K=25$, with bias	Training data II	-	98.36%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=25$, with bias	Training data II	-	98.37%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=10$	Training data II	-	98.01%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=20$	Training data II	-	98.23%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=40$	Training data II	-	98.22%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=60$	Training data II	-	98.25%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=80$	Training data II	-	98.22%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=10$, with bias	Training data II	-	98.23%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=20$, with bias	Training data II	-	98.46%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=40$, with bias	Training data II	-	98.54%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=60$, with bias	Training data II	-	98.60%
V	c	The Perceptron Criterion ($\eta(k)=\frac{1}{\sqrt{k}}$), $K=80$, with bias	Training data II	-	98.61%

Table 6: Part V classification results

In Part V, again, we do not know anything about the a priori knowledge. Although, all we know is the training data, we can train weight vectors $\vec{a}_i$ for linear discriminant functions $g_i(\vec{x})$ by minimizing a criterion function which in this case is the Perceptron criterion function. According to the results in the table above, the models achieve very promising results without knowing anything, i.e., this is actually a machine learning. Also, we can conclude that the problem is quite linearly separable since we only utilize the linear discriminant functions. The stopping criterion for all the Perceptron criterion above is stopping when the iteration $k$ is equal to the maximum number of iterations $K$.

We have tried a bunch of different $K$, and found out that $K=25$ is a good value to get good classification results. It is also interesting to see how the bias term $w_0$ would affect the performance of the models. With the same setting (i.e., the same $K$ and the same intial weights), we can see that the models with bias term $w_0$ yield better results in all cases.

For the question (c), we have the dynamic learning rate, which is a function of the current iteration, $\eta(k)$. This means the larger current iteration $k$ resulting in the smaller learning rate $\eta$, i.e., the learning rate decreases over time. In order to see how the dynamic learning rate has an effect to the performance of a model, we have varied $K \in \{10, 20, 40, 60, 80\}$. As can be seen from the table above, if we let a model trains further, the accuracy tend to increase. Nevertheless, one must be careful, if we train a model too much, it will eventually be overfit to the training data. There are many possible solutions to this issue, for instance, adding reguarlization, or applying cross validation.

---

All in all, this is actually a really wonderful MiniProject that wraps up pretty much everything about this course. Different supervised classification techniques with different knowledge and assumptions are employed in a series to see how each one works. The comparisons between them are also discussed. And this comes an end of this Introduction to Pattern Recognition and Machine Learning class.