{"id":22583,"date":"2021-02-01T11:43:00","date_gmt":"2021-02-01T11:43:00","guid":{"rendered":"https:\/\/www.experfy.com\/blog\/implement-expectation-maximization-algorithm-python\/"},"modified":"2023-09-05T11:45:51","modified_gmt":"2023-09-05T11:45:51","slug":"implement-expectation-maximization-algorithm-python","status":"publish","type":"post","link":"https:\/\/www.experfy.com\/blog\/ai-ml\/implement-expectation-maximization-algorithm-python\/","title":{"rendered":"Implement Expectation-Maximization Algorithm(EM) in Python from Scratch"},"content":{"rendered":"\t\t
\n\t\t\t\t\t\t
\n\t\t\t\t\t\t
\n\t\t\t\t\t
\n\t\t\t
\n\t\t\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

Unsupervised and Semi-supervised <\/strong>Gaussian Mixture Models (GMM)<\/p>\n

When companies launch a new product, they usually want to find out the target customers. If they have data on customers\u2019 purchasing history and shopping preferences, they can utilize it to predict what types of customers are more likely to purchase the new product. There are many models to solve this typical unsupervised learning problem and the Gaussian Mixture Model (GMM) is one of them.<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t

GMM and EM<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

GMMs are probabilistic models that assume all the data points are generated from a mixture of several Gaussian distributions with unknown parameters. They differ from k-means clustering in that GMMs incorporate information about the center(mean) and variability(variance) of each clusters and provide posterior probabilities.<\/p>\n

In the example mentioned earlier, we have 2 clusters: people who like the product and people who don\u2019t. If we know which cluster each customer belongs to (the labels), we can easily estimate the parameters(mean and variance) of the clusters, or if we know the parameters for both clusters, we can predict the labels. Unfortunately, we don\u2019t know either one. To solve this chicken and egg problem, the Expectation-Maximization Algorithm (EM) comes in handy.<\/p>\n

EM is an iterative algorithm to find the maximum likelihood when there are latent variables. The algorithm iterates between performing an expectation (E) step, which creates a heuristic of the posterior distribution and the log-likelihood using the current estimate for the parameters, and a maximization (M) step, which computes parameters by maximizing the expected log-likelihood from the E step. The parameter-estimates from M step are then used in the next E step. In the following sections, we will delve into the math behind EM, and implement it in Python from scratch.<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t

Mathematical Deduction<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

W define the known variables as x<\/em>, and the unknown label as y<\/em>. We make two assumptions: the prior distribution p(y)<\/em> is binomial and p(x|y) <\/em>in each cluster is a Gaussian .<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"variables\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

All parameters are randomly initialized. For simplicity, we use \u03b8 to represent all parameters in the following equations.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"\u03b8\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

At the expectation (E) step, we calculate the heuristics of the posteriors. We call them heuristics because they are calculated with guessed parameters \u03b8.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"E\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

At the maximization (M) step, we find the maximizers of the log-likelihood and use them to update \u03b8. Notice that the summation inside the logarithm in equation (3) makes the computational complexity NP-hard. To move the summation out of the logarithm, we use Jensen\u2019s inequality to find the evidence lower bound (ELBO) which is tight only when Q(y|x) = P(y|x). If you are interested in the math details from equation (3) to equation (5), this article<\/a> has decent explanation.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"Implement\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

Luckily, there are closed-form solutions for the maximizers in GMM.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"closed-form\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

We use these updated parameters in the next iteration of E step, get the new heuristics and run M-step. What the EM algorithm does is repeat these two steps until the average log-likelihood converges.<\/p>\n

Before jumping into the code, let\u2019s compare the above parameter solutions from EM to the direct parameter estimates when the labels are known. Did you find they are very similar? In fact, the only difference is that the EM solutions use the heuristics of posteriors Q<\/em> while the direct estimates use the true labels.<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"EM\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t

Python Implementation<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

There are many packages including scikit-learn that offer high-level APIs to train GMMs with EM. In this section, I will demonstrate how to implement the algorithm from scratch to solve both unsupervised and semi-supervised problems. The complete code can be find here<\/a>.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t

1. Unsupervised GMM<\/h3>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

Let\u2019s stick with the new product example. Using the known personal data, we have engineered 2 features x1, x2 <\/em>represented by a matrix x<\/em>, and our goal is to forecast whether each customer will like the product (y=1) or not (y=0).https:\/\/towardsdatascience.com\/media\/de0e60977c7b85502cad24a0abe46e8d<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t
\nhttps:\/\/gist.github.com\/VXU1230\/7333b0bf40256e1612d48af5080afea5\n<\/div><\/figure>\n

First we initialize all the unknown parameters.get_random_psd()<\/code> ensures the random initialization of the covariance matrices is positive semi-definite.<\/p>\n\n

https:\/\/gist.github.com\/VXU1230\/426921e27f0dfa6fa76aeaa95faec54f#file-initialize_random_params-py<\/p>\n

Then we pass the initialized parameters to e_step()<\/code>and calculate the heuristics Q(y=1|x) and Q(y=0|x)<\/em> for every data point as well as the average log-likelihoods which we will maximize in the M step.<\/p>\n\n

https:\/\/gist.github.com\/VXU1230\/dca3f37734759c278eaeb29a92020d80#file-e_step-py<\/p>\n

In m_step()<\/code> , the parameters are updated using the closed-form solutions in equation(7) ~ (11). Note that if there weren\u2019t closed-form solutions, we would need to solve the optimization problem using gradient ascent and find the parameter estimates.<\/p>\n\n

https:\/\/gist.github.com\/VXU1230\/3f6db97cce7a9c6969b83fa9c1c7d712#file-m_step-py<\/p>\n

Now we can repeat running the two steps until the average log-likelihood converges. rum_em()<\/code> returns the predicted labels, the posteriors and average log-likelihoods from all training steps.<\/p>\n\n

https:\/\/gist.github.com\/VXU1230\/c86fc58a1d3043a6fde6f48d6ccb1f9a#file-run_em-py<\/p>\n

Running the unsupervised model , we see the average log-likelihoods converged in over 30 steps.<\/p>\n\n

https:\/\/gist.github.com\/VXU1230\/ebf0771d2f9ef4635ad21eb1688dcd34#file-unsupervised-py<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"Implement\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t

2. Semi-supervised GMM<\/h3>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

In some cases, we have a small amount of labeled data. For example, we might know some customers\u2019 preferences from surveys. Considering the potential customer base is huge, the amount of labeled data we have is insufficient for full supervised learning, yet we can learn the initial parameters from the data in a semi-supervised way.<\/p>\n

We use the same unlabeled data as before, but we also have some labeled data this time.<\/p>

https:\/\/gist.github.com\/VXU1230\/9b6648f5202522c6ef2756e6ede74ab2#file-load_labeled_data-py<\/p>\n

In learn_params()<\/code> , we learn the initial parameters from the labeled data by implementing equation (12) ~(16). These learned parameters are used in the first E step.<\/p>

https:\/\/gist.github.com\/VXU1230\/77ca2a26190920d4724501453e84ba0c#file-learn_params-py<\/p>\n

Other than the initial parameters, everything else is the same so we can reuse the functions defined earlier. Let\u2019s train the model and plot the average log-likelihoods.<\/p>\n\n

https:\/\/gist.github.com\/VXU1230\/33b51d76fe736bd82d67a731ad50f441#file-semi-supervised-py<\/p>\n

This time the average log-likelihoods converged in 4 steps, much faster than unsupervised learning.<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"Implement\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

To verify our implementation, we compare our forecasts with forecasts from the scikit-learn API. To build the model in scikit-learn, we simply call the GaussianMixture API and fit the model with our unlabeled data. Don\u2019t forget to pass the learned parameters to the model so it has the same initialization as our semi-supervised implementation.GMM_sklearn()<\/code>returns the forecasts and posteriors from scikit-learn.<\/p>\n\n

https:\/\/gist.github.com\/VXU1230\/3e3dafedb33ed8b30689957ec356b969#file-compare_sklearn-py<\/p>\n

Comparing the results, we see that the learned parameters from both models are very close and 99.4% forecasts matched. In case you are curious, the minor difference is mostly caused by parameter regularization and numeric precision in matrix calculation.<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\"Implement\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

That\u2019s it! We just demystified the EM algorithm.<\/p>\n\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t\t

Conclusions<\/h2>\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t
\n\t\t\t\t
\n\t\t\t\t\t\t\t\t\t

In this article, we explored how to train Gaussian Mixture Models with the Expectation-Maximization Algorithm and implemented it in Python to solve unsupervised and semi-supervised learning problems. EM is a very useful method to find the maximum likelihood when the model depends on latent variables and therefore is frequently used in machine learning.<\/p>\n

I hope you had fun reading this article \ud83d\ude42<\/p>\t\t\t\t\t\t\t\t<\/div>\n\t\t\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/div>\n\t\t\t\t\t<\/div>\n\t\t<\/section>\n\t\t\t\t<\/div>\n\t\t","protected":false},"excerpt":{"rendered":"

This article explores how to train Gaussian Mixture Models with the Expectation-Maximization Algorithm and implemente it in Python to solve unsupervised and semi-supervised learning problems.<\/p>\n","protected":false},"author":1008,"featured_media":18537,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"content-type":"","footnotes":""},"categories":[183],"tags":[1283,92,1284,1269,1021],"ppma_author":[3744],"class_list":["post-22583","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-ai-ml","tag-gmm","tag-machine-learning","tag-semi-supervised-gmm","tag-semi-supervised-learning","tag-unsupervised-model"],"authors":[{"term_id":3744,"user_id":1008,"is_guest":0,"slug":"siwei-xu","display_name":"Siwei Xu","avatar_url":"https:\/\/www.experfy.com\/blog\/wp-content\/uploads\/2021\/05\/Siwei-Xu-150x150.jpeg","user_url":"https:\/\/www.linkedin.com\/company\/walmartglobaltech%20","last_name":"Xu","first_name":"Siwei","job_title":"","description":"Siwei Xu is Software Engineer at Walmart Labs."}],"_links":{"self":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts\/22583","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/users\/1008"}],"replies":[{"embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/comments?post=22583"}],"version-history":[{"count":4,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts\/22583\/revisions"}],"predecessor-version":[{"id":32334,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/posts\/22583\/revisions\/32334"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/media\/18537"}],"wp:attachment":[{"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/media?parent=22583"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/categories?post=22583"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/tags?post=22583"},{"taxonomy":"author","embeddable":true,"href":"https:\/\/www.experfy.com\/blog\/wp-json\/wp\/v2\/ppma_author?post=22583"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}