.. restindex page-title: Pyem crumb: Pyem link-title: Pyem encoding: utf-8 output-encoding: None tags: python,pyem,Expectation Maximization,EM,online EM,recursive EM file: basic_example1.py file: basic_example2.py file: basic_example3.py file: example1.png file: Bic_example.png /restindex .. Last Change: Mon Jul 02 09:00 PM 2007 J =================================================== PyEM, a python package for Gaussian mixture models =================================================== .. contents:: Tables of contents PyEM is a package which enables to create Gaussian Mixture Models (diagonal and full covariance matrices supported), to sample them, and to estimate them from data using Expectation Maximization algorithm. It can also draw confidence ellipsoides for multivariate models, and compute the Bayesian Information Criterion to assess the number of clusters in the data. In a near future, I hope to add so-called online EM (ie recursive EM) and variational Bayes implementation. PyEM is implemented in python, and uses the excellent numpy and scipy packages. Numpy is a python packages which gives python a fast multi-dimensional array capabilities (ala matlab and the likes); scipy leverages numpy to build common scientific features for signal processing, linear algebra, statistics, etc... Installation ============ .. _scipy: http://www.scipy.org Pyem depends on several packages to work: - numpy - matplotlib (if you wish to use the plotting facilities of pyem) Those packages are likely to be already installed in a typical numpy/scipy environment. Since september 2006, pyem is included in the sandbox of `scipy`_. The sandbox contains packages which are pending for approval in main scipy; that means it is not installed by default, and that you need to install scipy from sources. For the most up-to-date version of pyem, you need to download scipy from subversion, which contains the development branch of scipy. To install pyem, you just need to edit (or create if it does not exist) the file Lib/sandbox/enabled_packages.txt in scipy sources, and add one line with the name of the package (eg pyem). After, you just need to install scipy normally as explained `here `_. You can (and should) also test pyem installation using the following: .. raw:: html {+mycoloring} from scipy.sandbox import pyem pyem.test() {-mycoloring} basic usage ============ Once you are inside a python interpreter, you can import pyem using the follwing command: .. raw:: html {+mycoloring} from scipy.sandbox import pyem {-mycoloring} Creating, sampling and plotting a mixture ----------------------------------------- Importing pyem gives access to 3 classes: GM (for Gausssian Mixture), GMM (Gaussian Mixture Model) and EM (for Expectation Maximization). The first class GM can be used to create an artificial Gaussian Model, samples it, or plot it. The following example show how to create a 2 dimension Gaussian Model with 3 components, sample it and plot its confidence ellipsoids with matplotlib: .. raw:: html {mycolorize;input/softwares/pyem/basic_example1.py} .. raw:: latex \input{basic_example1.tex} which plots this figure: .. image:: example1.png :width: 500 :height: 400 There are basically two ways to create a GM instance: either an empty one (eg without mean, weights and covariances) by giving hyper parameters (dimension, number of clusters, type of covariance matrices) during instanciation, or giving all parameters using the class method GM.fromvalues. The first method is mostly useful as a container for learning. There are also methods to create random (but valid !) parameters for a Gaussian Mixture: this can be done by the function method GM.generate_params (a class method). Basic estimation of mixture parameters from data ------------------------------------------------ If you want to learn a Gaussian mixture from data with EM, you need to use two classes from pyem: GMM and EM. You first create a GMM object from a GM instance; then you can give the GMM object as a parameter to the EM class to compute iterations of EM; once the EM has finished the computation, the GM instance of GMM contains the computed parameters. .. raw:: html {mycolorize;input/softwares/pyem/basic_example2.py} .. raw:: latex \input{basic_example2.tex} GMM class do all the hard work for learning: it can compute the sufficient statistics given the current state of the model, and update its parameters from the sufficient statistics; the EM class uses a GMM instance to compute several iterations. The idea is that you can implements a different mixture model and uses the same EM class if you want (there are several optimized models for GMM, which are subclasses of GMM). Advanced topics =============== Bayesian Information Criterion and automatic clustering ------------------------------------------------------- The GMM class is also able to compute the bayesian information criterion (BIC), which can be used to assess the number of clusters into the data. It was first suggested by Schwarz (see bibliography), who gave a Bayesian argument for adopting the BIC. The BIC is derived from an approximation of the integrated likelihood of the model, based on regularity assumptions. The following code generates an artificial mixture of 4 clusters, runs EM with models of 1 to 6 clusters, and prints which number of clusters is the most likely from the BIC: .. raw:: html {mycolorize;input/softwares/pyem/basic_example3.py} .. raw:: latex \input{basic_example3.tex} which plots this figure: .. image:: Bic_example.png :width: 500 :height: 400 The above example also shows that you can control the plotting parameters by using returned handles from plot methods. This can be useful for complex drawing. Examples ========= Using EM for pdf estimation --------------------------- The following example uses the old faithful dataset and is available in the example directory. It models the joint distribution (d(t), w(t+1)), where d(t) is the duration time, and w(t+1) the waiting time for the next eruption. It selects the best model using the BIC. .. raw:: latex \input{pdfestimation.tex} .. figure:: pdfestimation.png :width: 500 :height: 400 isodensity curves for the old faithful data modeled by a 1, 2, 3 and 4 componenits model (up to bottom, left to right). Using EM for clustering ----------------------- TODO (this is fundamentally the same than pdf estimation, though) Using PyEM for supervised learning ---------------------------------- The following example shows how to do classification using discriminative analysis on the famous iris dataset. The idea is to model each class distribution P(data|class) by a mixture of Gaussian, and to classify the data with Maximum A Posteriori (eg takes the class which maximizes P(class | data)). Not surprisingly, the errors lie on the border between the two classes which are not linearly separable. .. raw:: latex \input{discriminant_analysis.tex} .. figure:: dclass.png :width: 600 :height: 400 learning data points on the left with isodensity curves as estimated by each mixture. On the right, the misclassified points are in red, green are correctly classified. Note on performances ==================== Pyem is implemented in python (100% of the code has a python implementation), but thanks to the Moore Law, it is reasonably fast so that it can be used for other problems than toys problem. On my computer (linux on bi xeon 3.2 Ghz, 2Gb RAM), running 10 iterations of EM algorithm on 100 000 samples of dimension 15, for a diagonal model with 30 components, takes around 1 minute and 15 seconds: this makes the implementation usable for moderately complex problems such as speaker recognition using MFCC. If this is too slow, there is a C implementation for Gaussian densities which can be enabled by commenting out one line in pyem/gmm_em.py, which should gives a speed up of a factor 2 at least; this has not been tested much, though, so beware. Also, increasing the number of components and/or dimension is much more expensive than increasing the number of samples: a 5 dimension model of 100 components will be much slower to estimate with 500 samples than a 5 dimension 10 components with 5000 samples. This is because loops on dimension/components are in python, whereas loops on samples are in C (through numpy). I don't think there is an easy fix to this problem. Full covariances will be slow, because you cannot avoid nested loop in python this case without insane amount of memory. A C implementation may be implemented, but this is not my top priority; most of the time, you should avoid full covariance models anyway. TODO ==== I believe the current API simple and powerful enough, except maybe for plotting (if you think otherwise, I would be happy to hear your suggestions). Now, I am considering adding some more functionalities to the toolbox: - add simple methods for regularization of covariance matrix (easy) - add bayes prior (using variational Bayes approximation) for overfitting and model selection problems (not trivial, but doable) - improve online EM Other things which are doable but which I don't intend to implement are: - add other models (mixtures of multinomial: easy, simple HMM: easy, other ?) - add bayes prior using MCMC (hard, use PyMCMC for sampling ?) Bibliography ============ TODO.