# EM

## Mixtures of Gaussians

given a training set $\left\{x^{(1)},...,x^{(n)}\right\}$ with no labels. each point $x^{(i)}$ has a latent variable $z^{(i)}$, the data is specified by a joint probability $p(x^{(i)}, z^{(i)}) = p(x^{(i)}| z^{(i)})p(z^{(i)})$

mixtures of gaussians assumes $z^{(i)}\sim Multinomial(\phi)$ (where $\phi_{j} \ge 0, \sum_{j=1}^{k}\phi_{j}=1$, $\phi_{j}$ gives $p(z^{(i)}=j)$), and $x^{(i)}|z^{(i)}=j \sim \mathcal{N}(\mu_{j}, \Sigma_{j})$. we let $k$ denote the number of values $z^{(i)}$'s can take on.

thus our model posits that each $x^{(i)}$ was generated by randomly choosing $z^{(i)}$ from $\{1,...,k\}$, and $x^{(i)}$ was drawn from one of $k$ gaussians depending on $z^{(i)}$.

the parameters of our model are $\phi, \mu, \Sigma$, to estimate them, we can write down the likelihood of our data:

$$
\begin{equation}
\begin{split}
l(\phi, \mu, \Sigma) =& \sum_{i=1}^{n}log\,p(x^{(i)};\phi,\mu,\Sigma)\\
=& \sum_{i=1}^{n}log\sum_{z^{(i)}=1}^{k}p(x^{(i)}|z^{(i)}; \mu,\Sigma)p(z^{(i)};\phi)
\end{split}
\end{equation}
$$

it is not possible to find the maximum likelihood estimate of the parameters in closed form.

on the other hand, if we knew what $z^{(i)}$'s were, the maximum likelihood problem would have been easy.

specifically, we could then write down the likelihood as:

$$l(\phi, \mu, \Sigma) = \sum_{i=1}^{n}log\,p(x^{(i)}|z^{(i)};\mu,\Sigma) + log\,p(z^{(i)}; \phi)$$

maximizing this with respect to $\phi, \mu, \Sigma$ gives the parameters:

$$\phi_{j} = \frac{1}{n}\sum_{i=1}^{n}1\{z^{(i)} = j\}$$

$$\mu_{j} = \frac{\sum_{i=1}^{n}1\{z^{(i)} = j\}x^{(i)}}{\sum_{i=1}^{n}1\{z^{(i)} = j\}}$$

$$\Sigma_{j} = \frac{\sum_{i=1}^{n}1\{z^{(i)} = j\}(x^{(i)} - \mu_{j})(x^{(i)} - \mu_{j})^{T}}{\sum_{i=1}^{n}1\{z^{(i)} = j\}}$$

it is identical to the gaussian discriminant analysis(GDA).

however, in our density estimation problems, $z^{(i)}$'s are not known, what can we do?

the EM algorithm is an iterative algorithm that has two main steps:

**1.the E-step, it tries to "guess" the values of the $z^{(i)}$'s.**

**2.the M-step, it updates the parameters of our model based on our guesses.**

since in the M-step we are pretending that the guesses in the first part were correct, the maximization becomes easy, here's the algorithm:


(E-step) for each $i, j$ set

$$w_{j}^{(i)} := p(z^{i}=j|x^{(i)}; \phi,\mu,\Sigma)$$

(M-step) update the parameters:

$$\phi_{j} := \frac{1}{n}\sum_{i=1}^{n}w_{j}^{(i)}$$

$$\mu_{j} := \frac{\sum_{i=1}^{n}w_{j}^{(i)}x^{(i)}}{\sum_{i=1}^{n}w_{j}^{(i)}}$$

$$\Sigma_{j} := \frac{\sum_{i=1}^{n}w_{j}^{(i)}(x^{(i)} - \mu_{j})(x^{(i)} - \mu_{j})^{T}}{\sum_{i=1}^{n}w_{j}^{(i)}}$$

the EM-algorithm likes the K-means-algorithm.

except that instead of "hard" cluster assignments $c(i)$, we have the "soft" assignments $w_{j}^{(i)}$.

from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs
import numpy as np

X1, y1 = make_blobs(n_samples=1000, centers=((4, -4), (0, 0)), random_state=42)
X1 = X1.dot(np.array([[0.374, 0.95], [0.732, 0.598]]))
X2, y2 = make_blobs(n_samples=250, centers=1, random_state=42)
X2 = X2 + [6, -8]
X = np.r_[X1, X2]
y = np.r_[y1, y2]

import matplotlib.pyplot as plt

def plot_clusters(X, y=None):
    plt.scatter(X[:, 0], X[:, 1], c=y, s=1)
    plt.xlabel("$x_1$", fontsize=14)
    plt.ylabel("$x_2$", fontsize=14, rotation=0)

plot_clusters(X)

gm = GaussianMixture(n_components=3, n_init=10, random_state=42)
gm.fit(X)

gm.weights_

gm.means_

gm.covariances_

gm.converged_, gm.n_iter_

gm.predict(X), gm.predict_proba(X)

X_new, y_new = gm.sample(7)

gm.score_samples(X_new)

gm.bic(X), gm.aic(X)

## The EM algorithm

here comes the general view of EM.

like the assumption before, we have the training set $\left\{x^{(1)},...,x^{(n)}\right\}$, each point $x^{(i)}$ has a latent variable $z^{(i)}$.

denote $\theta$ as the distribution parameters($\phi,\mu,\Sigma$ in mixtures of gaussian), then density of $x$ can be obtained by:

$$p(x;\theta) = \sum_{z}p(x,z;\theta)$$

the log-likelihood of the data:

$$
\begin{equation}
\begin{split}
l(\theta) =& \sum_{i=1}^{n}log\,p(x^{(i)}; \theta) \\
=& \sum_{i=1}^{n}log\sum_{z^{(i)}}p(x^{(i)}, z^{(i)}; \theta)
\end{split}
\end{equation}
$$

instead of solving the maximum-likelihood directly(usually very hard), the EM algorithm repeatedly construct a lower-bound on $l$(E-step), and then optimize that lower-bound(M-step).

we first consider optimizing the likelihood $log\,p(x)$ for a single example $x$, optimize $log\,p(x;\theta)$:

$$log\,p(x;\theta) = log\sum_{z}p(x,z;\theta)$$

let $Q$ be a distribution over the possible values of $z$, then:

$$
\begin{equation}
\begin{split}
log\,p(x;\theta) =& log\sum_{z}p(x,z;\theta)\\
=& log\sum_{z}Q(z)\frac{p(x,z;\theta)}{Q(z)}\\
\ge& \sum_{z}Q(z)log\frac{p(x,z;\theta)}{Q(z)}
\end{split}
\end{equation}
$$

the last step uses Jensen's inequality, for $(log\,x)'' = -1/x^{2} < 0$ so that strictly concave.

for any distribution $Q$, the above formula gives a lower-bound on $log\,p(x;\theta)$.

to make the bound tight for a particular value of $\theta$, we want the Jensen's inequality to hold with equality. it is sufficient the expecation take on "constant"-valued random variables, i.e:

$$\frac{p(x,z;\theta)}{Q(z)} = c$$

this leads to:

$$
\begin{equation}
\begin{split}
Q(z) =& \frac{p(x,z;\theta)}{\sum_{z}p(x,z;\theta)}\\
=& \frac{p(x,z;\theta)}{p(x;\theta)}\\
=& p(z|x;\theta)
\end{split}
\end{equation}
$$

thus, we simply set $Q$ to be the posterior distribution of $z$'s given $x$ and $\theta$.

for convenience, we call the bound as **evidence lower bound(ELBO)**:

$$ELBO(x;Q,\theta) = \sum_{z}Q(z)log\frac{p(x,z;\theta)}{Q(z)}$$

the estimation now can re-write as:

$$\forall Q,\theta,x,\quad log\,p(x;\theta) \ge ELBO(x;Q,\theta)$$

Intuitively, the EM algorithm alternatively update $Q$ and $\theta$ by:

a.setting $Q(z) = p(z|x;\theta)$ so that $ELBO(x;Q,\theta) = log\,p(x;\theta)$

b.maximizing $ELBO(x;Q,\theta)$ with respect to $\theta$.

now consider multiple examples $\left\{x^{(1)},...,x^{(n)}\right\}$, note the optimal choice of $Q$ is $p(z|x;\theta)$, and it depends on the particular example $x$.

therefore we introduce $n$ distributions $Q_{1},...,Q_{n}$, one for each example $x^{(i)}$, for each example:

$$log\,p(x^{(i)};\theta) \ge ELBO(x^{(i)};Q_{i},\theta) = \sum_{z^{(i)}}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$$

taking sum over all examples, we obtain a lower bound for the log-likelihood:

$$
\begin{equation}
\begin{split}
l(\theta) \ge& \sum_{i}ELBO(x^{(i)};Q_{i},\theta)\\
=& \sum_{i}\sum_{z^{(i)}}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}
\end{split}
\end{equation}
$$

the above inequality holds for **any** distributions $Q_{1},...,Q_{n}$. equality holds when:

$$Q_{i}(z^{(i)}) =  p(z^{(i)}|x^{(i)};\theta)$$

we now come to the definition of EM algorithm:

(E-step) for each $i$, set:

$$Q_{i}(z^{(i)}) :=  p(z^{(i)}|x^{(i)};\theta)$$

(M-step) for fixed $Q_{i}$'s, set:

$$
\begin{equation}
\begin{split}
\theta :=& \underset{\theta}{argmax}\sum_{i=1}^{n}ELBO(x^{(i)};Q_{i},\theta)\\
=& \underset{\theta}{argmax}\sum_{i=1}^{n}\sum_{z^{(i)}}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}
\end{split}
\end{equation}
$$

## Convergence of EM algorithm

suppose $\theta^{(t)}$ and $\theta^{(t+1)}$ are the parameters from two successive iterations of EM, we will now prove that $l(\theta^{(t)}) \le l(\theta^{(t+1)})$.

$$
\begin{equation}
\begin{split}
l(\theta^{(t+1)}) \ge& \sum_{i=1}^{n}ELBO(x^{(i)};Q_{i}^{(t)};\theta^{(t+1)})\quad (\mbox{by Jensen's inequality})\\
\ge& \sum_{i=1}^{n}ELBO(x^{(i)};Q_{i}^{(t)};\theta^{(t)})\quad (\mbox{by M-step}) \\
=& l(\theta^{(t)})\quad (\mbox{by E-step})
\end{split}
\end{equation}
$$

if we define:

$$ELBO(Q, \theta) = \sum_{i=1}^{n}ELBO(x^{(i)}; Q_{i}, \theta)$$

then we know $l(\theta) \ge ELBO(Q, \theta)$.

the EM algorithm can be viewed an alternating maximization algorithm on $ELBO(Q, \theta)$:

the E-step maximizes it with respect to $Q$.

the M-steo maximizes it with respect to $\theta$.

## Bayesian Gaussian Mixture Models

from sklearn.mixture import BayesianGaussianMixture

bgm = BayesianGaussianMixture(n_components=10, n_init=10)
bgm.fit(X)
np.round(bgm.weights_, 2)

## Excercises

### 1.clustering

from sklearn.datasets import fetch_olivetti_faces

olivetti = fetch_olivetti_faces()

print(olivetti.DESCR)

olivetti.target

from sklearn.model_selection import StratifiedShuffleSplit

strat_split = StratifiedShuffleSplit(n_splits=1, test_size=40, random_state=42)
train_valid_idx, test_idx = next(strat_split.split(olivetti.data, olivetti.target))
X_train_valid = olivetti.data[train_valid_idx]
y_train_valid = olivetti.target[train_valid_idx]
X_test = olivetti.data[test_idx]
y_test = olivetti.target[test_idx]

strat_split = StratifiedShuffleSplit(n_splits=1, test_size=80, random_state=43)
train_idx, valid_idx = next(strat_split.split(X_train_valid, y_train_valid))
X_train = X_train_valid[train_idx]
y_train = y_train_valid[train_idx]
X_valid = X_train_valid[valid_idx]
y_valid = y_train_valid[valid_idx]

print(X_train.shape, y_train.shape)
print(X_valid.shape, y_valid.shape)
print(X_test.shape, y_test.shape)

from sklearn.decomposition import PCA

pca = PCA(0.99)
X_train_pca = pca.fit_transform(X_train)
X_valid_pca = pca.transform(X_valid)
X_test_pca = pca.transform(X_test)

pca.n_components_

from sklearn.cluster import KMeans

k_range = range(5, 150, 5)
kmeans_per_k = []
for k in k_range:
    print("k={}".format(k))
    kmeans = KMeans(n_clusters=k, random_state=42).fit(X_train_pca)
    kmeans_per_k.append(kmeans)

from sklearn.metrics import silhouette_score

silhouette_scores = [silhouette_score(X_train_pca, model.labels_)
                     for model in kmeans_per_k]
best_index = np.argmax(silhouette_scores)
best_k = k_range[best_index]
best_score = silhouette_scores[best_index]

plt.figure(figsize=(8, 3))
plt.plot(k_range, silhouette_scores, "bo-")
plt.xlabel("$k$", fontsize=14)
plt.ylabel("Silhouette score", fontsize=14)
plt.plot(best_k, best_score, "rs")
plt.show()

best_k

inertias = [model.inertia_ for model in kmeans_per_k]
best_inertia = inertias[best_index]

plt.figure(figsize=(8, 3.5))
plt.plot(k_range, inertias, "bo-")
plt.xlabel("$k$", fontsize=14)
plt.ylabel("Inertia", fontsize=14)
plt.plot(best_k, best_inertia, "rs")
plt.show()

best_model = kmeans_per_k[best_index]

def plot_faces(faces, labels, n_cols=5):
    faces = faces.reshape(-1, 64, 64)
    n_rows = (len(faces) - 1) // n_cols + 1
    plt.figure(figsize=(n_cols, n_rows * 1.1))
    for index, (face, label) in enumerate(zip(faces, labels)):
        plt.subplot(n_rows, n_cols, index + 1)
        plt.imshow(face, cmap="gray")
        plt.axis("off")
        plt.title(label)
    plt.show()

for cluster_id in np.unique(best_model.labels_):
    print("Cluster", cluster_id)
    in_cluster = best_model.labels_==cluster_id
    faces = X_train[in_cluster]
    labels = y_train[in_cluster]
    plot_faces(faces, labels)

### 2.clustering + classification

from sklearn.ensemble import RandomForestClassifier

clf = RandomForestClassifier(n_estimators=150, random_state=42)
clf.fit(X_train_pca, y_train)
clf.score(X_valid_pca, y_valid)

X_train_reduced = best_model.transform(X_train_pca)
X_valid_reduced = best_model.transform(X_valid_pca)
X_test_reduced = best_model.transform(X_test_pca)

clf = RandomForestClassifier(n_estimators=150, random_state=42)
clf.fit(X_train_reduced, y_train)
    
clf.score(X_valid_reduced, y_valid)

from sklearn.pipeline import Pipeline

for n_clusters in k_range:
    pipeline = Pipeline([
        ("kmeans", KMeans(n_clusters=n_clusters, random_state=42)),
        ("forest_clf", RandomForestClassifier(n_estimators=150, random_state=42))
    ])
    pipeline.fit(X_train_pca, y_train)
    print(n_clusters, pipeline.score(X_valid_pca, y_valid))

X_train_extended = np.c_[X_train_pca, X_train_reduced]
X_valid_extended = np.c_[X_valid_pca, X_valid_reduced]
X_test_extended = np.c_[X_test_pca, X_test_reduced]

clf = RandomForestClassifier(n_estimators=150, random_state=42)
clf.fit(X_train_extended, y_train)
clf.score(X_valid_extended, y_valid)

### 3.Gaussian Mixture Model

from sklearn.mixture import GaussianMixture

gm = GaussianMixture(n_components=40, random_state=42)
y_pred = gm.fit_predict(X_train_pca)

n_gen_faces = 20
gen_faces_reduced, y_gen_faces = gm.sample(n_samples=n_gen_faces)
gen_faces = pca.inverse_transform(gen_faces_reduced)

plot_faces(gen_faces, y_gen_faces)

n_rotated = 4
rotated = np.transpose(X_train[:n_rotated].reshape(-1, 64, 64), axes=[0, 2, 1])
rotated = rotated.reshape(-1, 64*64)
y_rotated = y_train[:n_rotated]

n_flipped = 3
flipped = X_train[:n_flipped].reshape(-1, 64, 64)[:, ::-1]
flipped = flipped.reshape(-1, 64*64)
y_flipped = y_train[:n_flipped]

n_darkened = 3
darkened = X_train[:n_darkened].copy()
darkened[:, 1:-1] *= 0.3
y_darkened = y_train[:n_darkened]

X_bad_faces = np.r_[rotated, flipped, darkened]
y_bad = np.concatenate([y_rotated, y_flipped, y_darkened])

plot_faces(X_bad_faces, y_bad)

X_bad_faces_pca = pca.transform(X_bad_faces)

gm.score_samples(X_bad_faces_pca)

gm.score_samples(X_train_pca[:10])

### 4.Using Dimensionality Reduction Techniques for Anomaly Detection

X_train_pca

def reconstruction_errors(pca, X):
    X_pca = pca.transform(X)
    X_reconstructed = pca.inverse_transform(X_pca)
    mse = np.square(X_reconstructed - X).mean(axis=-1)
    return mse

reconstruction_errors(pca, X_train).mean()

reconstruction_errors(pca, X_bad_faces).mean()

plot_faces(X_bad_faces, y_bad)

X_bad_faces_reconstructed = pca.inverse_transform(X_bad_faces_pca)
plot_faces(X_bad_faces_reconstructed, y_bad)

