# Bernoulli Mixture

### Bernoulli Mixture following Bishop’s Pattern Recognition and Machine Learning Section 9.3.3.

The code below gives a basic implementation of the Bernoulli Mixture fit using the EM algorithm. The essential equations are 9.57 and 9.58 in Bishop. We apply it to the fashion-mnist dataset which has 10 classes, but we look for 20 classes and pick up variations among the purses, etc. The resulting figure is basically a version of Figure 9.10 in Bishop. Different runs of this code produce different results!

# %load bernoulli.py
import pandas as pd
import numpy as np
from numpy import log
from scipy.special import logsumexp
import matplotlib
import matplotlib.pyplot as plt
from itertools import product

class BernoulliMixture:

def __init__(self, n_clusters=2, n_iter=100, tolerance=1e-8, alpha1=1e-6, alpha2=1e-6):
'''sets things up'''
self._n_clusters = n_clusters
self._n_iter = n_iter
self._tolerance = tolerance
self._alpha1 = alpha1
self._alpha2 = alpha2
self.Theta = 1/self._n_clusters*np.ones(self._n_clusters)

def _P(self, x, Mu, Theta):
'''computes the log of the conditional probability of the latent variable given the data and Mu, Theta'''
ll = np.dot(x,log(Mu))+np.dot(1-x,log(1-Mu))
Z = (log(Theta)+ ll - logsumexp(ll+log(Theta), axis=1, keepdims=True))
return Z

def fit(self, data):
'''carries out the EM iteration'''
self._n_samples, self._n_features = data.shape
self.Mu = np.random.uniform(.25, .75, size=self._n_clusters*self._n_features).reshape(self._n_features, self._n_clusters)
N = self._n_samples

for i in range(self._n_iter):
V = np.exp(self._P(data, self.Mu, self.Theta))
W = V/V.sum(axis=1,keepdims=True)
R = np.dot(data.transpose(), W)
Q = np.sum(W, axis=0, keepdims=True)
self._old = self.Theta
self.Mu, self.Theta = (R + self._alpha1)/(Q + self._n_features*self._alpha1), (Q + self._alpha2)/(N + self._n_clusters * self._alpha2)
if np.allclose(self._old, self.Theta):
return

def predict_proba(self, data):
'''computes the conditional probability giving cluster membership'''
return np.exp(self._P(data, self.Mu, self.Theta))

def generate(self):
'''generates data from the distribution'''
x = np.random.multinomial(1,self.Theta.ravel())
z = np.dot(np.random.binomial(1, self.Mu),self.Theta.transpose())
return z

if __name__ == "__main__":
data = data_df.iloc[:,1:].values
data = data // 128
M = BernoulliMixture(n_clusters=20)
M.fit(data)
fig, ax = plt.subplots(5,4)
fig.set_size_inches(12,12)
for i,j in product(range(5), range(4)):
if 4*i+j<20:
ax[i,j].imshow(M.Mu[:, 4*i+j].reshape(28,28))
else:
ax[i,j].axis('off')
fig.savefig('trial.png') Updated: