# Gaussian Mixture Implementation and Optical Recognition of Handwritten Digits Data Set

Trying to implement Gaussian Mixture model implementation in python using the Optical Recognition of Handwritten Digits Data Set which consists of 10 training folds each of size $$\left[100x64\right]$$, and 10 training labels each of size $$\left[100x1\right]$$. The dataset also has a testing dataset and label set of size $$\left[110x64\right]$$, and $$\left[110x1\right]$$. There are only two classes 5 and 6. I get the following error on the class conditional density method:

ValueError: operands could not be broadcast together with shapes (64,64) (100,64)


According to Bishop (Pattern Recognition And Machine Learning 2006). I first have to estimate optimal parameters through MLE. So for every fold I'm estimating the optimal parameters, but then I don't know how to calculate posterior probabilities given 10 different folds and 10 labels. $$p\left(x\mid C_{k}\right)=\frac{1}{\left(2\pi\right)^{\frac{D}{2}}}\cdot\frac{1}{\left|\sum\right|^{\frac{1}{2}}}\cdot exp\left\{ -\frac{1}{2}\left(x-\mu_{k}\right)^{T}\cdot\sum^{-1}\cdot\left(x-\mu_{k}\right)\right\}$$ $$pi^{*}=\frac{N_{1}}{N_{1}+N_{2}}$$

$$mu_{1}^{*}=\frac{1}{N_{1}}\sum_{n=1}^{N}t_{n}x_{n}$$, $$\mu_{2}=\frac{1}{N_{2}}\sum_{n=1}^{N}t_{n}x_{n}$$

$$cov^{*}=-\frac{N}{2}ln\left|\Sigma\right|-\frac{N}{2}Tr\left\{ \Sigma^{-1}S\right\}$$

$$S^{*}=\frac{N_{1}}{N}s_{1}+\frac{N_{2}}{N}s_{2}$$

$$s_{1}^{*}=\frac{1}{N_{1}}\sum_{n\in C_{1}}\left(x_{n}-\mu_{1}\right)\left(x_{n}-\mu_{1}\right)^{T}$$

$$s_{2}^{*}=\frac{1}{N_{2}}\sum_{n\in C_{2}}\left(x_{n}-\mu_{2}\right)\left(x_{n}-\mu_{2}\right)^{T}$$

Then, I have to calculate posterior probabilities with $$p\left(C_{1}\mid x\right)=\sigma\left(w^{T}x+w_{0}\right)$$

$$a_{k}\left(x\right)=w_{k}^{T}x+w_{k0}$$

$$w_{k}=\Sigma^{-1}\mu_{k}$$

$$w_{k0}=-\frac{1}{2}\mu_{k}^{T}\Sigma^{-1}\mu_{k}+lnp\left(C_{k}\right)$$ I don't even know if my approach is correct. I've searched on either GitHub and medium for similar examples. Any help or guidance would be appreciated. My code implementation:

#----------------------------------------------------------
# Assigning parameters
#----------------------------------------------------------
#--------------------------------------------------------
# Initial guess
#----------------------------------------------------------
# Total Folds
folds = 10
#Define the number of classes
num_classes = 2
n, m = total_trainf[0].shape
means = np.zeros((num_classes , n))
phi = np.zeros((num_classes , 1))
shared_cov_matrix = np.cov(np.transpose(np.concatenate(total_trainf, axis=0)),bias=True)
#posterior_prob_est = pp(w0, w, np.concatenate(total_trainf, axis=1))
c= 1
f= 0
acc_arry = []
# calculate the maximum likelihood of each observation xi
likelihood = []
means_array = np.zeros([num_classes, total_trainf[0].shape[1]])
cond_class_prob = np.zeros([num_classes, 1])
# Expectation step
means = np.zeros([folds, total_trainf[0].shape[1]])
while f < folds:
uclasses = np.unique(total_trainl[f])
mu_arr, pi, sigma = get_params(total_trainf[f], total_trainl[f], means_array, phi, uclasses)
#means[f, :] = mu_arr
cov = get_cov(total_trainl[f], shared_cov_matrix, sigma)
class_prob = ccd(total_trainf[f], mu_arr, pi, shared_cov_matrix )
f+=1
means_array[c-1] = np.mean(means, axis=0)

def get_params(data, label, means, pi, class_type):
mean_array = np.zeros([1, data.shape[1]])
num_classes = len(class_type)
sigma = 0
col = 0
for i in range(num_classes):
ind = np.flatnonzero(label == class_type[i])
pi[i] = len(ind)/label.shape[0]
means[i] = np.mean(data[ind] , axis = 0)
sigma += np.cov(data[ind].T)*(len(ind) - 1)

sigma = sigma/label.shape[0]
return means, pi, sigma

#Covariance per fold
def get_cov(data, cov, S):
res = np.log(np.linalg.det(cov))+np.trace(np.linalg.inv(cov)*S)
return -(data.shape[0]/2)*res

#Class conditional density
def ccd(data, means: np.array, pi: np.array, cov_matrix: np.array):
inv_cov = np.linalg.inv(cov_matrix)
mu_1 = means[0]
mu_2 = means[1]
pi_1 = pi[0]
pi_2 = pi[1]
W = inv_cov*(mu_1 - mu_2)
ft = -0.5*(mu_1.T*inv_cov*mu_1)
st = 0.5*(mu_2.T*inv_cov*mu_2)
tt = np.log(pi_1/pi_2)
W_0 = ft + st + tt
return pp(W, W_0, data)
$$$$
`