{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 1. Introduction\n\n
\n

1.1 Objectives

\n
\n\nIn this notebook we'll be building a **Gaussian Mixture Model** (GMM), **Poisson Mixture Model** (PMM) and a **Hybrid Mixture Model** (HMM) from scratch. We'll use the **Expectation-Maximisation** (EM) algorithm to learn the parameters of these models.\n\n**Note:** There is very **little information about PMMs** online and it seems they haven't been studied much before. As a result, a lot of the equations in this notebook have been derived by myself; some rigorously and others using only my intuition of EM. I am also quite confident that there will be more efficient implementations of these models. If you have any ideas for improvements please share your feedback. \n\nMy other [notebook on GMMs](https://www.kaggle.com/code/samuelcortinhas/gaussian-mixture-model-gmm-from-the-ground-up) goes into more detail of how the EM algorithm works, so feel free to read that one before this one.\n\n**TL;DR:** The data for this competition isn't modelled well by PMM or HMM but the code can still be used for other problems.\n\n
\n

1.2 Libraries

\n
","metadata":{}},{"cell_type":"code","source":"# Core\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nsns.set(style='darkgrid', font_scale=1.4)\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom matplotlib import cm\nfrom itertools import combinations\nimport math\nimport statistics\nfrom scipy import stats\nfrom scipy.stats import pearsonr\nfrom scipy.stats import shapiro\nfrom scipy.stats import chi2\nfrom scipy.stats import poisson\nfrom scipy.stats import multivariate_normal\nfrom scipy.special import factorial\nfrom scipy.stats import poisson\nfrom scipy.stats import norm\nimport time\nfrom datetime import datetime\nimport matplotlib.dates as mdates\nimport plotly.express as px\nfrom termcolor import colored\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\n# Sklearn\nfrom sklearn.decomposition import PCA\nfrom sklearn.manifold import TSNE\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\nfrom sklearn.cluster import KMeans\nfrom sklearn.model_selection import train_test_split, StratifiedKFold, GridSearchCV, TimeSeriesSplit\nfrom sklearn.preprocessing import StandardScaler, RobustScaler, PowerTransformer, OneHotEncoder, LabelEncoder\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.compose import make_column_transformer\nfrom sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, accuracy_score\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.linear_model import LinearRegression, Ridge\nfrom sklearn.mixture import GaussianMixture, BayesianGaussianMixture\n\n# UMAP\nimport umap\nimport umap.plot","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-17T10:14:00.654557Z","iopub.execute_input":"2022-07-17T10:14:00.655664Z","iopub.status.idle":"2022-07-17T10:14:36.571819Z","shell.execute_reply.started":"2022-07-17T10:14:00.655530Z","shell.execute_reply":"2022-07-17T10:14:36.570585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Gaussian Mixture Model\n\n
\n

2.1 Gaussian Mixture

\n
\n\nA **Gaussian Mixture** is simply a combination (or mixture) of Gaussian distributions. In **d dimensions**, the model is a weighted sum of **multivariate** normal distributions:\n\n$$\nf_{\\text{GMM}} (\\textbf{x}) = \\sum_{j=1}^{k} \\pi_j f_{\\mathcal{N}({\\boldsymbol \\mu}_j, {\\boldsymbol \\Sigma}_{j})} (\\textbf{x})\n$$\n\nwhere\n\n* $\\textbf{x} = (x_1, \\ldots, x_d)$ is a vector of length $d$\n* $f_{\\mathcal{N}({\\boldsymbol \\mu}_j, {\\boldsymbol \\Sigma}_{j})}$ is the density of a **multivariate** normal distribution with mean vector ${\\boldsymbol \\mu}$ and covariance matrix ${\\boldsymbol \\Sigma}$\n* $\\pi = (\\pi_1, \\ldots, \\pi_k)$ are the weights subject to\n\n$$\n0 \\leq \\pi_j \\leq 1, \\quad \\sum_{j=1}^{k} \\pi_j = 1.\n$$","metadata":{}},{"cell_type":"code","source":"# Define grid\nx = np.linspace(0,20, num=100).astype(int)\ny = np.linspace(0,20, num=100).astype(int)\nx, y = np.meshgrid(x, y)\n\n# Calculate pdf over mesh\nz1 = multivariate_normal.pdf(np.dstack((x, y)), [6,5], [[4,1],[-1,4]])\nz2 = multivariate_normal.pdf(np.dstack((x, y)), [14,14], [[3,1],[1,3]])\nz3 = multivariate_normal.pdf(np.dstack((x, y)), [16,7], [[10,1],[-1,10]])\nz = z1 + z2 + z3\n\n# Plot 3D plot with gaussian mixture\nfig = plt.figure(figsize=(6,6))\nax = fig.add_subplot(111, projection='3d')\nax.plot_surface(x,y,z, cmap=cm.jet)\nax.set_title('2D Gaussian mixture plotted in 3D')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-17T10:14:36.575551Z","iopub.execute_input":"2022-07-17T10:14:36.577366Z","iopub.status.idle":"2022-07-17T10:14:37.481591Z","shell.execute_reply.started":"2022-07-17T10:14:36.577299Z","shell.execute_reply":"2022-07-17T10:14:37.480380Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

2.2 EM algorithm for the GMM

\n
\n\n\n## Initialisation\n\nFor each cluster $j$, choose the **mean** ${\\boldsymbol \\mu}_j$ to be a **random data point** and the **covariance** ${\\boldsymbol \\Sigma}_j$ to be the **covariance of the whole dataset** $X$. The **weights** $\\pi = (\\pi_1, \\ldots, \\pi_k)$ are initially **uniform**.\n\n## E-step\n\nWe update the **responsibilities** (i.e. posterior probabilities) using Bayes' formula, where $r_{ij}$ is the probability that the $i$-th data point belongs to the $j$-th mixture\n\n$$\nr_{ij} = \\mathbb{P}(C_j | \\textbf{x}_i) = \\frac{\\mathbb{P}(\\textbf{x}_i|C_j) \\mathbb{P}(C_j)}{\\sum_{t=1}^{k} \\mathbb{P}(\\textbf{x}_i|C_t) \\mathbb{P}(C_t)} = \\frac{\\mathbb{P}(\\textbf{x}_i|C_j) \\pi_j}{\\sum_{t=1}^{k} \\mathbb{P}(\\textbf{x}_i|C_t) \\pi_t}\n$$\n\nwhere the **likelihoods** are Gaussian densities\n\n$$\n\\mathbb{P}(\\textbf{x}_i | C_j) = \\frac{1}{\\sqrt{(2 \\pi)^d \\det({\\boldsymbol \\Sigma}_j)}} \\exp \\left(-\\frac{(\\textbf{x}_i-{\\boldsymbol \\mu}_j)^T {\\boldsymbol \\Sigma}_j^{-1} (\\textbf{x}_i-{\\boldsymbol \\mu}_j)}{2} \\right)\n$$\n\n## M-step\n\nWe estimate the **mean** and **covariances** of each class as follows\n\n$$\n{\\boldsymbol \\mu}_j = \\frac{\\sum_{i=1}^n r_{ij} \\textbf{x}_{i}}{n_j}, \\qquad {\\boldsymbol \\Sigma}_j = \\frac{1}{n_j} \\sum_{i=1}^{n} r_{ij} (\\textbf{x}_i - {\\boldsymbol \\mu}_j) (\\textbf{x}_i - {\\boldsymbol \\mu}_j)^{T}\n$$\n\nwhere $n_j = \\sum_{i=1}^{n} r_{ij}$ is defined as the **total responsibility** of the j-th mixture component.\n\nThe **mixture weights** are updated via\n\n$$\n\\pi_j = \\frac{n_j}{n}.\n$$\n\n
\n

2.3 GMM implementation

\n
","metadata":{}},{"cell_type":"code","source":"class GMM:\n def __init__(self, k, max_iter=100, random_state = 0):\n self.k = k\n self.max_iter = max_iter\n self.random_state = random_state\n\n def initialise(self, X):\n self.shape = X.shape\n self.n, self.d = self.shape\n \n self.pi = np.full(shape=self.k, fill_value=1/self.k)\n self.responsibilities = np.full(shape=self.shape, fill_value=1/self.k)\n \n np.random.seed(self.random_state)\n random_row = np.random.randint(low=0, high=self.n, size=self.k)\n self.mu = [X[row_index,:] for row_index in random_row]\n self.sigma = [np.cov(X.T) for _ in range(self.k)]\n\n def E_step(self, X):\n # E-Step: update the responsibilities by holding mu and sigma constant\n self.responsibilities = self.predict_proba(X)\n \n def M_step(self, X):\n # M-Step: update pi, mu and sigma by holding responsibilities constant\n self.pi = self.responsibilities.mean(axis=0)\n for j in range(self.k):\n r_column = self.responsibilities[:,j]\n total_responsibility = r_column.sum()\n self.mu[j] = (X * r_column[:, np.newaxis]).sum(axis=0)/total_responsibility\n self.sigma[j] = np.cov(X.T, aweights=(r_column/total_responsibility).flatten(), bias=True)\n\n def fit(self, X):\n self.initialise(X)\n \n for iteration in range(self.max_iter):\n self.E_step(X)\n self.M_step(X)\n \n def predict_proba(self, X):\n likelihood = np.zeros((self.n, self.k))\n for j in range(self.k):\n distribution = multivariate_normal(mean=self.mu[j], cov=self.sigma[j])\n likelihood[:,j] = distribution.pdf(X)\n \n numerator = likelihood * self.pi\n denominator = numerator.sum(axis=1)[:, np.newaxis]\n responsibilities = numerator / denominator\n return responsibilities\n \n def predict(self, X):\n responsibilities = self.predict_proba(X)\n return np.argmax(responsibilities, axis=1)\n \n def fit_predict(self, X):\n self.fit(X)\n predictions = self.predict(X)\n return predictions","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T20:57:17.651772Z","iopub.execute_input":"2022-07-16T20:57:17.652264Z","iopub.status.idle":"2022-07-16T20:57:17.675018Z","shell.execute_reply.started":"2022-07-16T20:57:17.652231Z","shell.execute_reply":"2022-07-16T20:57:17.674180Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Poisson Mixture Model\n\n
\n

3.1 Poisson Mixture

\n
\n\nA **Poisson Mixture** is a mixture of Poisson distributions. In **d dimensions**, the model is a weighted sum of **multivariate** Poisson distributions:\n\n$$\nf_{\\text{PMM}} (\\textbf{x}) = \\sum_{j=1}^{k} \\pi_j f_{\\text{PoM}({\\boldsymbol \\Lambda}_j)} (\\textbf{x})\n$$\n\nwhere the d-dimensional **multivariate Poisson distribution** corresponds to the **discrete** [distribution](https://reference.wolfram.com/language/ref/MultivariatePoissonDistribution.html#:~:text=The%20multivariate%20Poisson%20distribution%20is,and%20covariance%20of%20the%20distribution.)\n\n$$\n(X_0 + X_1, X_0 + X_2, \\ldots, X_0 + X_d)\n$$\n\nwith\n\n$$\nX_i \\sim \\text{Po}(\\lambda_j).\n$$\n\n
\n\nRecall that $X \\sim \\text{Po}(\\lambda)$ means that\n\n$$\n\\mathbb{P}(X=k) = \\frac{e^{-\\lambda} \\lambda^{k}}{k!}\n$$\n\n
\n\nWe start by deriving the **probability mass function** of the multivariate poisson distribution. This will be our **likelihood** in the EM algorithm. \n\nNotice that the event $(X_0 + X_1, X_0 + X_2, \\ldots, X_0 + X_d) = (x_1, x_2, \\ldots, x_d)$ is the **disjoint union** of the events $(X_0, X_1, \\ldots, X_d) = (i, x_1 - i, x_2 -i, \\ldots, x_d - i)$.\n\nThis means we can **sum** over the probabilities that are well defined, which happens when $i \\leq \\min(x_1, x_2, \\ldots, x_d)$, otherwise the values would be negative. This gives\n\n$$\n\\begin{align*}\nf_{\\text{PoM}({\\boldsymbol \\Lambda})} (\\textbf{x}) &= \\sum_{i=0}^{\\min(x_1, \\ldots, x_d)} \\mathbb{P}(X_0 = i) \\mathbb{P}(X_1 = x_1-i) \\ldots \\mathbb{P}(X_d = x_d -i) \\newline\n& = \\sum_{i=0}^{\\min(x_1, \\ldots, x_d)} \\left(\\frac{e^{-\\lambda_0} \\lambda_0^{i}}{i!} \\right) \\left(\\frac{e^{-\\lambda_1} \\lambda_1^{(x_1 - i)}}{(x_1 - i)!} \\right) \\ldots \\left(\\frac{e^{-\\lambda_d} \\lambda_d^{(x_d - i)}}{(x_d - i)!} \\right) \\newline\n&= e^{-(\\lambda_0 + \\ldots + \\lambda_d)} \\sum_{i=0}^{\\min(x_1, \\ldots, x_d)} \\frac{\\lambda_0^{i}}{i!} \\frac{\\lambda_1^{(x_1 - i)}}{(x_1 - i)!} \\ldots \\frac{\\lambda_d^{(x_d - i)}}{(x_d - i)!}\n\\end{align*}\n$$\n\n(I think that's the simplest form.)\n","metadata":{}},{"cell_type":"code","source":"# Define grid\nx = np.linspace(0,25, num=26).astype(int)\ny = np.linspace(0,25, num=26).astype(int)\nx, y = np.meshgrid(x, y)\nz = 0\n\n# Calculate pmf\ndef multivariate_poisson(X, Y, mu0, mu1, mu2):\n exponent = np.exp(-(mu0+mu1+mu2))\n out = np.zeros(X.shape)\n for p in range(len(X)):\n for q in range(len(X)):\n x = X[p,q]\n y = Y[p,q]\n minimum = np.min((x,y))\n term=0\n for i in range(minimum+1):\n powers = np.array([i, x-i, y-i])\n term += np.prod(np.power([mu0, mu1, mu2],powers) / factorial(powers, exact=True))\n out[p,q] = exponent * term\n return out\n\n# Evaluate pmf over grid\nz1 = multivariate_poisson(x, y, 1, 2, 3)\nz2 = multivariate_poisson(x, y, 10, 6, 3)\nz = z1 + z2\n\n# Plot 3D plot with gaussian mixture\nfig = plt.figure(figsize=(6,6))\nax = fig.add_subplot(111, projection='3d')\nax.plot_surface(x,y,z, cmap=cm.jet)\nax.set_title('2D Poisson mixture plotted in 3D')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T20:57:17.678695Z","iopub.execute_input":"2022-07-16T20:57:17.679392Z","iopub.status.idle":"2022-07-16T20:57:18.742854Z","shell.execute_reply.started":"2022-07-16T20:57:17.679359Z","shell.execute_reply":"2022-07-16T20:57:18.741831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

3.2 EM algorithm for the PMM

\n
\n\n\n## Initialisation\n\nFor each cluster $j$, we choose the **vector** $(\\lambda_1, \\ldots, \\lambda_d)$ to be a **random data point** and $\\lambda_0$ to be a random interger between 0 and k, which together make up the multivariate **rate** parameter $\\Lambda$. The **weights** $\\pi = (\\pi_1, \\ldots, \\pi_k)$ are initially **uniform**.\n\n## E-step\n\nWe update the **responsibilities** (i.e. posterior probabilities) using Bayes' formula, where $r_{ij}$ is the probability that the $i$-th data point belongs to the $j$-th mixture\n\n$$\nr_{ij} = \\mathbb{P}(C_j | \\textbf{x}_i) = \\frac{\\mathbb{P}(\\textbf{x}_i|C_j) \\mathbb{P}(C_j)}{\\sum_{t=1}^{k} \\mathbb{P}(\\textbf{x}_i|C_t) \\mathbb{P}(C_t)} = \\frac{\\mathbb{P}(\\textbf{x}_i|C_j) \\pi_j}{\\sum_{t=1}^{k} \\mathbb{P}(\\textbf{x}_i|C_t) \\pi_t}\n$$\n\nwhere the **likelihoods** are multivariate Poisson densities\n\n$$\n\\mathbb{P}(C_j | \\textbf{x}_i) = e^{-(\\lambda_0 + \\ldots + \\lambda_d)} \\sum_{i=0}^{\\min(x_1, \\ldots, x_d)} \\frac{\\lambda_0^{i}}{i!} \\frac{\\lambda_1^{(x_1 - i)}}{(x_1 - i)!} \\ldots \\frac{\\lambda_d^{(x_d - i)}}{(x_d - i)!}\n$$\n\n## M-step\n\nWe estimate the **mean** and **covariances** of each class as follows\n\n$$\n{\\boldsymbol \\mu}_j = \\frac{\\sum_{i=1}^n r_{ij} \\textbf{x}_{i}}{n_j}, \\qquad {\\boldsymbol \\Sigma}_j = \\frac{1}{n_j} \\sum_{i=1}^{n} r_{ij} (\\textbf{x}_i - {\\boldsymbol \\mu}_j) (\\textbf{x}_i - {\\boldsymbol \\mu}_j)^{T}\n$$\n\nwhere $n_j = \\sum_{i=1}^{n} r_{ij}$ is defined as the **total responsibility** of the j-th mixture component.\n\n
\n\nIt is known [theory](https://reference.wolfram.com/language/ref/MultivariatePoissonDistribution.html#:~:text=The%20multivariate%20Poisson%20distribution%20is,and%20covariance%20of%20the%20distribution.) that the **covariance matrix** of the multivariate Poisson distribution has the following form\n\n$$\n\\left(\\begin{array}{cccc}\n(\\lambda_0 + \\lambda_1) & \\lambda_0 & \\ldots & \\lambda_0 \\\\\n\\lambda_0 & (\\lambda_0 + \\lambda_2) & {} & \\vdots \\\\\n\\vdots & {} & \\ddots & \\lambda_0 \\\\\n\\lambda_0 & \\ldots & \\lambda_0 & (\\lambda_0 + \\lambda_d)\n\\end{array}\\right)\n$$\n\nIn practice, the non-diagonal elements are not all the same because of the noise in the dataset. What I propose, is to **average** all the **non-diagonal** elements to calculate $\\lambda_0$, and then to update the Poisson rates via ${\\boldsymbol \\mu} - \\lambda_0 {\\boldsymbol 1}$.\n\n
\n\nThe **mixture weights** are updated via\n\n$$\n\\pi_j = \\frac{n_j}{n}.\n$$\n\n
\n

3.3 PMM implementation

\n
","metadata":{}},{"cell_type":"code","source":"class PMM:\n def __init__(self, k, max_iter=10, random_state=0):\n self.k = k\n self.max_iter = max_iter\n self.random_state = random_state\n\n def initialise(self, X):\n self.shape = X.shape\n self.n, self.d = self.shape\n \n self.pi = np.full(shape=self.k, fill_value=1/self.k)\n self.responsibilities = np.full(shape=self.shape, fill_value=1/self.k)\n \n np.random.seed(self.random_state)\n lambda0 = np.random.choice(range(self.k), size=self.k, replace = False)\n random_row = np.random.choice(range(self.n), size=self.k, replace = False)\n self.lambd = np.c_[lambda0, np.array([X[row_index,:] for row_index in random_row])] # shape (k,d+1)\n \n def E_step(self, X):\n # E-Step: update the responsibilities by holding lambda constant\n self.responsibilities = self.predict_proba(X)\n \n def M_step(self, X):\n # M-Step: update pi and lambda by holding responsibilities constant\n self.pi = self.responsibilities.mean(axis=0)\n for j in range(self.k):\n r_column = self.responsibilities[:,j]\n r_column[r_column<0]=0\n total_responsibility = r_column.sum()\n mu = (X * r_column[:, np.newaxis]).sum(axis=0)/total_responsibility\n sigma = np.cov(X.T, aweights=(r_column/total_responsibility).flatten(), bias=True)\n lambda0 = (sigma.sum()-np.diag(sigma).sum())/(self.d**2-self.d) # average non-diagonal elements\n self.lambd[j,:] = np.concatenate(([lambda0], mu-lambda0))\n \n def fit(self, X):\n self.initialise(X)\n \n for iteration in range(self.max_iter):\n self.E_step(X)\n self.M_step(X)\n \n def multivariate_poisson(self, X, cluster):\n exponent = np.exp(-self.lambd[cluster,:].sum())\n probs=[]\n for p in range(self.n):\n x = X[p,:]\n term=0\n for i in range(x.min()+1):\n powers = np.concatenate(([i],[x[q]-i for q in range(self.d)]))\n term += np.prod(np.power(self.lambd[cluster,:],powers)/factorial(powers, exact=True))\n probs.append(exponent * term)\n return probs\n \n def predict_proba(self, X):\n likelihood = np.zeros((self.n, self.k))\n for j in range(self.k):\n likelihood[:,j] = self.multivariate_poisson(X,j)\n \n numerator = likelihood * self.pi\n denominator = numerator.sum(axis=1)[:, np.newaxis]\n responsibilities = numerator / denominator\n return responsibilities\n \n def predict(self, X):\n responsibilities = self.predict_proba(X)\n return np.argmax(responsibilities, axis=1)\n \n def fit_predict(self, X):\n self.fit(X)\n predictions = self.predict(X)\n return predictions","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T20:57:18.745396Z","iopub.execute_input":"2022-07-16T20:57:18.746204Z","iopub.status.idle":"2022-07-16T20:57:18.776998Z","shell.execute_reply.started":"2022-07-16T20:57:18.746162Z","shell.execute_reply":"2022-07-16T20:57:18.775033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Hybrid Mixture Model\n\n
\n

4.1 Hybrid Mixture

\n
\n\nWe define a **Hybrid Mixture** to be the mixture of hybrid distributions. For the purposes of the [July 2022 TPS competition](https://www.kaggle.com/competitions/tabular-playground-series-jul-2022), we define the relevant hybrid distributions to be the **concatenation** of a **multivariate Poisson distribution** with a **multivariate Gaussian distribution**. In particular,\n\n$$\n(X, Y)\n$$\n\nwhere\n\n$$\nX = (X_0+X_1, X_0+X_2, \\ldots, X_0+X_{d_1})\n$$\n\n$$\nX_i \\sim \\, \\text{Po}(\\lambda_i), \\quad i=0, \\ldots, d_1\n$$\n\nand\n\n$$\nY = (Y_1, Y_2, \\ldots, Y_{d_2})\n$$\n\n$$\nY_i \\sim \\, \\mathcal{N}({\\boldsymbol \\mu}_i,{\\boldsymbol \\Sigma}_i), \\quad i=1, \\ldots, d_2.\n$$\n\n
\n\nFor $f_{\\text{Hyb}} \\sim (X,Y)$, our **Hybrid Mixture Model** is \n\n$$\nf_{\\text{HMM}} (\\textbf{x}, \\textbf{y}) = \\sum_{j=1}^{k} \\pi_j f_{\\text{Hyb}({\\boldsymbol \\Lambda}_j, {\\boldsymbol \\mu}_j, {\\boldsymbol \\Sigma}_j)} (\\textbf{x}, \\textbf{y}) = \\sum_{j=1}^{k} \\pi_j (f_{\\text{PoM}({\\boldsymbol \\Lambda}_j)} (\\textbf{x}), f_{\\mathcal{N}({\\boldsymbol \\mu}_j, {\\boldsymbol \\Sigma}_j)} (\\textbf{y})).\n$$","metadata":{}},{"cell_type":"markdown","source":"
\n

4.2 EM algorithm for the HMM

\n
\n\n\n## Initialisation\n\nFor each cluster $j$, we initialise the **Poisson parameters** like in a PMM and the **Gaussian parameters** like in a GMM (see above). The **weights** $\\pi = (\\pi_1, \\ldots, \\pi_k)$ are initially **uniform**.\n\n## E-step\n\nWe update the **responsibilities** (i.e. posterior probabilities) using Bayes' formula, where $r_{ij}$ is the probability that the $i$-th data point belongs to the $j$-th mixture\n\n$$\nr_{ij} = \\mathbb{P}(C_j | (\\textbf{x}, \\textbf{y})_i) = \\frac{\\mathbb{P}((\\textbf{x}, \\textbf{y})_i|C_j) \\mathbb{P}(C_j)}{\\sum_{t=1}^{k} \\mathbb{P}((\\textbf{x}, \\textbf{y})_i|C_t) \\mathbb{P}(C_t)} = \\frac{\\mathbb{P}((\\textbf{x}, \\textbf{y})_i|C_j) \\pi_j}{\\sum_{t=1}^{k} \\mathbb{P}((\\textbf{x}, \\textbf{y})_i|C_t) \\pi_t}\n$$\n\nwhere the **likelihoods** are the **product** of the multivariate Poisson and multivariate Gaussian mass/density functions \n\n$$\n\\mathbb{P}(C_j | (\\textbf{x}, \\textbf{y})_i) = f_{\\text{PoM}({\\boldsymbol \\Lambda}_j)} (\\textbf{x}_i) f_{\\mathcal{N}({\\boldsymbol \\mu}_j, {\\boldsymbol \\Sigma}_{j})} (\\textbf{y}_i)\n$$\n\nThe key idea is that the Poisson parameters and the Gaussian parameters **share** the same responsibilities.\n\n## M-step\n\nFor the **Gaussian parameters**\n\n$$\n{\\boldsymbol \\mu}_j = \\frac{\\sum_{i=1}^n r_{ij} \\textbf{y}_{i}}{n_j}, \\qquad {\\boldsymbol \\Sigma}_j = \\frac{1}{n_j} \\sum_{i=1}^{n} r_{ij} (\\textbf{y}_i - {\\boldsymbol \\mu}_j) (\\textbf{y}_i - {\\boldsymbol \\mu}_j)^{T}\n$$\n\nwhere $n_j = \\sum_{i=1}^{n} r_{ij}$ is defined as the **total responsibility** of the j-th mixture component.\n\n
\n\nFor the **Poisson parameters**\n\n$$\n\\hat{{\\boldsymbol \\mu}}_j = \\frac{\\sum_{i=1}^n r_{ij} \\textbf{x}_{i}}{n_j}, \\qquad \\hat{{\\boldsymbol \\Sigma}}_j = \\frac{1}{n_j} \\sum_{i=1}^{n} r_{ij} (\\textbf{x}_i - \\hat{{\\boldsymbol \\mu}}_j) (\\textbf{x}_i - \\hat{{\\boldsymbol \\mu}}_j)^{T}\n$$\n\nwhere $n_j = \\sum_{i=1}^{n} r_{ij}$ is defined as the **total responsibility** of the j-th mixture component.\n\nThen **average** all the **non-diagonal** elements of $\\hat{{\\boldsymbol \\Sigma}}_j$ to calculate $\\lambda_0$, and then to update the Poisson rates via $\\hat{{\\boldsymbol \\mu}} - \\lambda_0 {\\boldsymbol 1}$.\n\n
\n\nThe **mixture weights** are updated via\n\n$$\n\\pi_j = \\frac{n_j}{n}.\n$$\n\n
\n

4.3 HMM implementation

\n
","metadata":{}},{"cell_type":"code","source":"class HMM:\n def __init__(self, k, d1, max_iter=10, random_state=0, verbose=True):\n self.k = k\n self.d1 = d1\n self.max_iter = max_iter\n self.random_state = random_state\n self.verbose = True\n\n def initialise(self, X):\n self.shape = X.shape\n self.n, self.d = self.shape # d1+d2=d\n \n self.pi = np.full(shape=self.k, fill_value=1/self.k)\n self.responsibilities = np.full(shape=self.shape, fill_value=1/self.k)\n \n np.random.seed(self.random_state)\n lambda0 = np.random.choice(range(self.k), size=self.k, replace = False)\n random_row = np.random.choice(range(self.n), size=self.k, replace = False)\n self.lambd = np.c_[lambda0, np.array([X[row_index,:self.d1] for row_index in random_row])] # shape (k,d1+1)\n \n self.mu = [X[row_index,self.d1:] for row_index in random_row] # shape (k,d2)\n self.sigma = [np.cov(X[:,self.d1:].T) for _ in range(self.k)] # shape (k,d2,d2)\n \n def E_step(self, X):\n # E-Step: update the responsibilities by holding lambda, mu, sigma constant\n self.responsibilities = self.predict_proba(X)\n \n def M_step(self, X):\n # M-Step: update pi, lambda, mu and sigma by holding responsibilities constant\n self.pi = self.responsibilities.mean(axis=0)\n for j in range(self.k):\n # Shared responsibilities\n r_column = self.responsibilities[:,j]\n r_column[r_column<0]=0\n total_responsibility = r_column.sum()\n \n # Poisson part\n mu_poisson = (X[:,:self.d1] * r_column[:, np.newaxis]).sum(axis=0)/total_responsibility\n sigma_poisson = np.cov(X[:,:self.d1].T, aweights=(r_column/total_responsibility).flatten(), bias=True)\n lambda0 = (sigma_poisson.sum()-np.diag(sigma_poisson).sum())/(self.d1**2-self.d1) # average non-diagonal elements\n self.lambd[j,:] = np.concatenate(([lambda0], mu_poisson-lambda0))\n \n # Gaussian part\n self.mu[j] = (X[:,self.d1:] * r_column[:, np.newaxis]).sum(axis=0)/total_responsibility\n self.sigma[j] = np.cov(X[:,self.d1:].T, aweights=(r_column/total_responsibility).flatten(), bias=True)\n\n def fit(self, X):\n self.initialise(X)\n \n for iteration in range(self.max_iter):\n if self.verbose == True:\n print('Iteration: ',iteration)\n self.E_step(X)\n self.M_step(X)\n \n def multivariate_poisson(self, X, cluster):\n exponent = np.exp(-self.lambd[cluster,:].sum())\n probs=[]\n for p in range(self.n):\n x = X[p,:]\n term=0\n for i in range(x.min()+1):\n powers = np.concatenate(([i],[x[q]-i for q in range(self.d1)]))\n term += np.prod(np.power(self.lambd[cluster,:],powers)/factorial(powers, exact=True))\n probs.append(exponent * term)\n return probs\n \n def predict_proba(self, X):\n likelihood = np.zeros((self.n, self.k))\n for j in range(self.k):\n gaussian = multivariate_normal(mean=self.mu[j], cov=self.sigma[j])\n likelihood[:,j] = self.multivariate_poisson(X[:,:self.d1].astype(int),j) * gaussian.pdf(X[:,self.d1:])\n \n numerator = likelihood * self.pi\n denominator = numerator.sum(axis=1)[:, np.newaxis]\n responsibilities = numerator / denominator\n return responsibilities\n \n def predict(self, X):\n responsibilities = self.predict_proba(X)\n return np.argmax(responsibilities, axis=1)\n \n def fit_predict(self, X):\n self.fit(X)\n predictions = self.predict(X)\n return predictions","metadata":{"execution":{"iopub.status.busy":"2022-07-16T20:57:18.779003Z","iopub.execute_input":"2022-07-16T20:57:18.779454Z","iopub.status.idle":"2022-07-16T20:57:18.849563Z","shell.execute_reply.started":"2022-07-16T20:57:18.779412Z","shell.execute_reply":"2022-07-16T20:57:18.848370Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

4.4 Train HMM

\n
","metadata":{}},{"cell_type":"code","source":"%%time\n\n# Load and preprocess data\ndata=pd.read_csv('../input/tabular-playground-series-jul-2022/data.csv', index_col='id')\ndrop_feats = [f'f_0{i}' for i in range(7)]\ndrop_feats = drop_feats + [f'f_{i}' for i in range(14,22)]\nX = data.drop(drop_feats, axis=1).values\n\n# Hybrid Mixture Model\nhmm = HMM(k=7, d1=7, max_iter=25, random_state=42)\ny = hmm.fit_predict(X)","metadata":{"execution":{"iopub.status.busy":"2022-07-16T20:57:18.851086Z","iopub.execute_input":"2022-07-16T20:57:18.851527Z","iopub.status.idle":"2022-07-16T21:05:12.938460Z","shell.execute_reply.started":"2022-07-16T20:57:18.851486Z","shell.execute_reply":"2022-07-16T21:05:12.937371Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot label distribution\nplt.figure(figsize=(10,4))\nsns.countplot(y)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

4.5 Evaluate predictions

\n
","metadata":{}},{"cell_type":"markdown","source":"**Continuous features**","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/code/ambrosm/tpsjul22-gaussian-mixture-cluster-analysis\nfig, axs = plt.subplots(2, 4, figsize=(20, 7))\naxs = axs.ravel()\nfloat_columns = ['f_22','f_23','f_24','f_25','f_26','f_27','f_28']\nfor ax, f in zip(axs, float_columns):\n for i in range(7):\n h, edges = np.histogram(data[f][y == i], bins=np.linspace(-5, 5, 26))\n ax.plot((edges[:-1] + edges[1:]) / 2, h, label=f\"Cluster {i}\", lw=3)\n ax.set_title(f)\naxs[-1].axis('off')\nplt.suptitle('Histograms of continuous features by cluster', y=1.02, fontsize=28)\nfig.tight_layout(h_pad=1.0, w_pad=0.5)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T21:05:12.939555Z","iopub.execute_input":"2022-07-16T21:05:12.939831Z","iopub.status.idle":"2022-07-16T21:05:14.165848Z","shell.execute_reply.started":"2022-07-16T21:05:12.939801Z","shell.execute_reply":"2022-07-16T21:05:14.164842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Normal Q-Q plots**","metadata":{}},{"cell_type":"code","source":"# Normal Q-Q plots\nfor i, col in enumerate(['f_22','f_23','f_24','f_25','f_26','f_27','f_28']):\n fig = plt.figure(figsize=(30,4))\n for j in range(7):\n clus = data[col][y==j]\n ax = plt.subplot(1, 7, j+1)\n stats.probplot(clus, dist='norm', plot=plt)\n \n # Aesthetics\n ax.get_lines()[0].set_markersize(6.0)\n ax.get_lines()[1].set_linewidth(3.0)\n ax.get_lines()[0].set_markerfacecolor(f'C{j}')\n ax.set_ylabel('')\n ax.set_xlabel('')\n if j==0:\n ax.set_ylabel('Observed values')\n if i==6:\n ax.set_xlabel('Theoretical quantiles')\n ax.set_xticklabels([])\n ax.set_yticklabels([])\n plt.title(f'cluster={j}')\n fig.tight_layout(h_pad=1.0, w_pad=0.5)\n plt.suptitle(col, y=1.05, fontsize=28)\n plt.show()\n print('')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T21:05:14.167407Z","iopub.execute_input":"2022-07-16T21:05:14.167709Z","iopub.status.idle":"2022-07-16T21:05:21.277072Z","shell.execute_reply.started":"2022-07-16T21:05:14.167681Z","shell.execute_reply":"2022-07-16T21:05:21.275691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Discrete features**","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/code/ambrosm/tpsjul22-gaussian-mixture-cluster-analysis\nprop_cycle = plt.rcParams['axes.prop_cycle']\n\nfig, axs = plt.subplots(2, 4, figsize=(20, 7))\naxs = axs.ravel()\nint_columns = [col for col in data.columns if data[col].dtype == 'int']\nfor ax, f in zip(axs, int_columns):\n for i in range(7):\n uv, uc = np.unique(data[f][y == i], return_counts=True)\n ax.plot(uv, uc, alpha=1, color=prop_cycle.by_key()['color'][i % 10], lw=3)\n ax.set_title(f)\n #ax.legend()\naxs[-1].axis('off')\nplt.suptitle('Histograms of discrete features by cluster', y=1.02, fontsize=28)\nfig.tight_layout(h_pad=1.0, w_pad=0.5)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T21:05:21.280139Z","iopub.execute_input":"2022-07-16T21:05:21.280449Z","iopub.status.idle":"2022-07-16T21:05:22.411903Z","shell.execute_reply.started":"2022-07-16T21:05:21.280424Z","shell.execute_reply":"2022-07-16T21:05:22.409878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Poisson Q-Q plots**","metadata":{}},{"cell_type":"code","source":"# Poisson Q-Q plots\nfor i, col in enumerate(['f_07','f_08','f_09','f_10','f_11','f_12','f_13']):\n fig = plt.figure(figsize=(30,4))\n for j in range(7):\n clus = data[col][y==j]\n mu = clus.mean()\n ax = plt.subplot(1, 7, j+1)\n stats.probplot(clus, dist='poisson', sparams=(mu,), plot=plt)\n \n # Aesthetics\n ax.get_lines()[0].set_markersize(6.0)\n ax.get_lines()[1].set_linewidth(3.0)\n ax.get_lines()[0].set_markerfacecolor(f'C{j}')\n ax.set_ylabel('')\n ax.set_xlabel('')\n if j==0:\n ax.set_ylabel('Observed values')\n if i==6:\n ax.set_xlabel('Theoretical quantiles')\n ax.set_xticklabels([])\n ax.set_yticklabels([])\n plt.title(f'cluster={j}')\n fig.tight_layout(h_pad=1.0, w_pad=0.5)\n plt.suptitle(col, y=1.05, fontsize=28)\n plt.show()\n print('')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T21:05:22.413465Z","iopub.execute_input":"2022-07-16T21:05:22.413813Z","iopub.status.idle":"2022-07-16T21:05:32.369390Z","shell.execute_reply.started":"2022-07-16T21:05:22.413782Z","shell.execute_reply":"2022-07-16T21:05:32.368549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**PCA**","metadata":{}},{"cell_type":"code","source":"%%time\n\n# PCA\npca = PCA(n_components=3)\ncomponents = pca.fit_transform(X)\n\n# 3D scatterplot\nfig = px.scatter_3d(\n components, x=0, y=1, z=2, color=y, size=0.1*np.ones(len(X)), opacity = 1,\n title='PCA plot in 3D',\n labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'},\n width=650, height=500\n)\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-16T21:05:32.370600Z","iopub.execute_input":"2022-07-16T21:05:32.371125Z","iopub.status.idle":"2022-07-16T21:05:34.188737Z","shell.execute_reply.started":"2022-07-16T21:05:32.371096Z","shell.execute_reply":"2022-07-16T21:05:34.187877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n

4.6 Submit predictions

\n
","metadata":{}},{"cell_type":"code","source":"sub = pd.read_csv('../input/tabular-playground-series-jul-2022/sample_submission.csv')\nsub['Predicted'] = y\nsub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2022-07-16T21:05:34.189795Z","iopub.execute_input":"2022-07-16T21:05:34.190216Z","iopub.status.idle":"2022-07-16T21:05:34.384297Z","shell.execute_reply.started":"2022-07-16T21:05:34.190189Z","shell.execute_reply":"2022-07-16T21:05:34.382744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5. Conclusion\n\nUnfortunately, the Poisson and Hybrid Mixture Models **don't perform well** on the leaderboard. As pointed out in this [discussion post](https://www.kaggle.com/competitions/tabular-playground-series-jul-2022/discussion/337563), the discrete features often have **negative covariances** which cannot be modelled properly using the multivariate Poisson distribution we defined above. The data in general has a weird distribution and it appears to be better modelled by scalling and then using a GMM. \n\nEven though the PMM and HMM didn't perform well on this dataset, I am still happy to have been able to derive and implement them from scratch since there are no other implementations available online. Hopefully it can be useful to someone else in the future.\n\n# 6. References\n\nThese will be resources I used when creating the PMM. For GMM resources see my other notebook linked right at the top. \n\n* [Multivariate Poisson Distribution](https://reference.wolfram.com/language/ref/MultivariatePoissonDistribution.html#:~:text=The%20multivariate%20Poisson%20distribution%20is,and%20covariance%20of%20the%20distribution.) by Wolfram.\n* [Poisson Distribution](https://en.wikipedia.org/wiki/Poisson_distribution) by Wikipedia.\n* [Bivariate Poisson Distribution](https://stats.stackexchange.com/questions/108705/deriving-the-bivariate-poisson-distribution) by whuber. \n* [Univariate Poisson Mixture Model](https://www.cs.helsinki.fi/u/bmmalone/probabilistic-models-spring-2014/PoissonMixtureModels.pdf) by Brandon Malone. ","metadata":{}}]}