مدل آمیخته گوسی در خوشه بندی – پیاده‌سازی در پایتون

1import numpy as np # linear algebra
2import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
3from sklearn.cluster import KMeans 
4from sklearn.feature_extraction.text import TfidfVectorizer
5from sklearn.decomposition import PCA
6from sklearn.preprocessing import normalize
7from sklearn.metrics import pairwise_distances
8
9import nltk
10import string
11
12import matplotlib.pyplot as plt
13%matplotlib inline
14plt.style.use('fivethirtyeight')
15
16
17# email module has some useful functions
18import os, sys, email,re
19
20df = pd.read_csv('../input/emails.csv',nrows = 35000)
21
22sklearn_pca = PCA(n_components = 2)
23Y_sklearn = sklearn_pca.fit_transform(tf_idf_array)
24kmeans = KMeans(n_clusters=3, max_iter=600, algorithm = 'auto')
25fitted = kmeans.fit(Y_sklearn)
26prediction = kmeans.predict(Y_sklearn)

1plt.figure(figsize = (10,8))
2from scipy.spatial.distance import cdist
3def plot_kmeans(kmeans, X, n_clusters=3, rseed=0, ax=None):
4    labels = kmeans.fit_predict(X)
5
6    # plot the input data
7    ax = ax or plt.gca()
8    ax.axis('equal')
9    ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)
10
11    # plot the representation of the KMeans model
12    centers = kmeans.cluster_centers_
13    radii = [cdist(X[labels == i], [center]).max()
14             for i, center in enumerate(centers)]
15    for c, r in zip(centers, radii):
16        ax.add_patch(plt.Circle(c, r, fc='#CCCCCC', lw=3, alpha=0.5, zorder=1))
17        
18plot_kmeans(kmeans, Y_sklearn)

1class GMM:
2    """ Gaussian Mixture Model
3    
4    Parameters
5    -----------
6        k: int , number of gaussian distributions
7        
8        seed: int, will be randomly set if None
9        
10        max_iter: int, number of iterations to run algorithm, default: 200
11        
12    Attributes
13    -----------
14       centroids: array, k, number_features
15       
16       cluster_labels: label for each data point
17       
18    """
19    def __init__(self, C, n_runs):
20        self.C = C # number of Guassians/clusters
21        self.n_runs = n_runs
22        
23    
24    def get_params(self):
25        return (self.mu, self.pi, self.sigma)
26    
27    
28        
29    def calculate_mean_covariance(self, X, prediction):
30        """Calculate means and covariance of different
31            clusters from k-means prediction
32        
33        Parameters:
34        ------------
35        prediction: cluster labels from k-means
36        
37        X: N*d numpy array data points 
38        
39        Returns:
40        -------------
41        intial_means: for E-step of EM algorithm
42        
43        intial_cov: for E-step of EM algorithm
44        
45        """
46        d = X.shape[1]
47        labels = np.unique(prediction)
48        self.initial_means = np.zeros((self.C, d))
49        self.initial_cov = np.zeros((self.C, d, d))
50        self.initial_pi = np.zeros(self.C)
51        
52        counter=0
53        for label in labels:
54            ids = np.where(prediction == label) # returns indices
55            self.initial_pi[counter] = len(ids[0]) / X.shape[0]
56            self.initial_means[counter,:] = np.mean(X[ids], axis = 0)
57            de_meaned = X[ids] - self.initial_means[counter,:]
58            Nk = X[ids].shape[0] # number of data points in current gaussian
59            self.initial_cov[counter,:, :] = np.dot(self.initial_pi[counter] * de_meaned.T, de_meaned) / Nk
60            counter+=1
61        assert np.sum(self.initial_pi) == 1    
62            
63        return (self.initial_means, self.initial_cov, self.initial_pi)

1class GMM:
2    """ Gaussian Mixture Model
3    
4    Parameters
5    -----------
6        k: int , number of gaussian distributions
7        
8        seed: int, will be randomly set if None
9        
10        max_iter: int, number of iterations to run algorithm, default: 200
11        
12    Attributes
13    -----------
14       centroids: array, k, number_features
15       
16       cluster_labels: label for each data point
17       
18    """
19    def __init__(self, C, n_runs):
20        self.C = C # number of Guassians/clusters
21        self.n_runs = n_runs
22        
23    
24    def get_params(self):
25        return (self.mu, self.pi, self.sigma)
26    
27    
28        
29    def calculate_mean_covariance(self, X, prediction):
30        """Calculate means and covariance of different
31            clusters from k-means prediction
32        
33        Parameters:
34        ------------
35        prediction: cluster labels from k-means
36        
37        X: N*d numpy array data points 
38        
39        Returns:
40        -------------
41        intial_means: for E-step of EM algorithm
42        
43        intial_cov: for E-step of EM algorithm
44        
45        """
46        d = X.shape[1]
47        labels = np.unique(prediction)
48        self.initial_means = np.zeros((self.C, d))
49        self.initial_cov = np.zeros((self.C, d, d))
50        self.initial_pi = np.zeros(self.C)
51        
52        counter=0
53        for label in labels:
54            ids = np.where(prediction == label) # returns indices
55            self.initial_pi[counter] = len(ids[0]) / X.shape[0]
56            self.initial_means[counter,:] = np.mean(X[ids], axis = 0)
57            de_meaned = X[ids] - self.initial_means[counter,:]
58            Nk = X[ids].shape[0] # number of data points in current gaussian
59            self.initial_cov[counter,:, :] = np.dot(self.initial_pi[counter] * de_meaned.T, de_meaned) / Nk
60            counter+=1
61        assert np.sum(self.initial_pi) == 1    
62            
63        return (self.initial_means, self.initial_cov, self.initial_pi)

1from scipy.stats import multivariate_normal as mvn

1      def _e_step(self, X, pi, mu, sigma):
2        """Performs E-step on GMM model
3
4        Parameters:
5        ------------
6        X: (N x d), data points, m: no of features
7        pi: (C), weights of mixture components
8        mu: (C x d), mixture component means
9        sigma: (C x d x d), mixture component covariance matrices
10
11        Returns:
12        ----------
13        gamma: (N x C), probabilities of clusters for objects
14        """
15        N = X.shape[0] 
16        self.gamma = np.zeros((N, self.C))
17
18        const_c = np.zeros(self.C)
19        
20        
21        self.mu = self.mu if self._initial_means is None else self._initial_means
22        self.pi = self.pi if self._initial_pi is None else self._initial_pi
23        self.sigma = self.sigma if self._initial_cov is None else self._initial_cov
24
25        for c in range(self.C):
26            # Posterior Distribution using Bayes Rule
27            self.gamma[:,c] = self.pi[c] * mvn.pdf(X, self.mu[c,:], self.sigma[c])
28
29        # normalize across columns to make a valid probability
30        gamma_norm = np.sum(self.gamma, axis=1)[:,np.newaxis]
31        self.gamma /= gamma_norm
32
33        return self.gamma

1      def _m_step(self, X, gamma):
2        """Performs M-step of the GMM
3        We need to update our priors, our means
4        and our covariance matrix.
5
6        Parameters:
7        -----------
8        X: (N x d), data 
9        gamma: (N x C), posterior distribution of lower bound 
10
11        Returns:
12        ---------
13        pi: (C)
14        mu: (C x d)
15        sigma: (C x d x d)
16        """
17        N = X.shape[0] # number of objects
18        C = self.gamma.shape[1] # number of clusters
19        d = X.shape[1] # dimension of each object
20
21        # responsibilities for each gaussian
22        self.pi = np.mean(self.gamma, axis = 0)
23
24        self.mu = np.dot(self.gamma.T, X) / np.sum(self.gamma, axis = 0)[:,np.newaxis]
25
26        for c in range(C):
27            x = X - self.mu[c, :] # (N x d)
28            
29            gamma_diag = np.diag(self.gamma[:,c])
30            x_mu = np.matrix(x)
31            gamma_diag = np.matrix(gamma_diag)
32
33            sigma_c = x.T * gamma_diag * x
34            self.sigma[c,:,:]=(sigma_c) / np.sum(self.gamma, axis = 0)[:,np.newaxis][c]
35
36        return self.pi, self.mu, self.sigma
37    
38    
39    def _compute_loss_function(self, X, pi, mu, sigma):
40        """Computes lower bound loss function
41        
42        Parameters:
43        -----------
44        X: (N x d), data 
45        
46        Returns:
47        ---------
48        pi: (C)
49        mu: (C x d)
50        sigma: (C x d x d)
51        """
52        N = X.shape[0]
53        C = self.gamma.shape[1]
54        self.loss = np.zeros((N, C))
55
56        for c in range(C):
57            dist = mvn(self.mu[c], self.sigma[c],allow_singular=True)
58            self.loss[:,c] = self.gamma[:,c] * (np.log(self.pi[c]+0.00001)+dist.logpdf(X)-np.log(self.gamma[:,c]+0.000001))
59        self.loss = np.sum(self.loss)
60        return self.loss

1      def fit(self, X):
2        """Compute the E-step and M-step and
3            Calculates the lowerbound
4        
5        Parameters:
6        -----------
7        X: (N x d), data 
8        
9        Returns:
10        ----------
11        instance of GMM
12        
13        """
14        
15        d = X.shape[1]
16        self.mu, self.sigma, self.pi =  self._initialise_parameters(X)
17        
18        try:
19            for run in range(self.n_runs):  
20                self.gamma  = self._e_step(X, self.mu, self.pi, self.sigma)
21                self.pi, self.mu, self.sigma = self._m_step(X, self.gamma)
22                loss = self._compute_loss_function(X, self.pi, self.mu, self.sigma)
23                
24                if run % 10 == 0:
25                    print("Iteration: %d Loss: %0.6f" %(run, loss))
26
27        
28        except Exception as e:
29            print(e)
30            
31        
32        return self

1     def predict(self, X):
2        """Returns predicted labels using Bayes Rule to
3        Calculate the posterior distribution
4        
5        Parameters:
6        -------------
7        X: N*d numpy array
8        
9        Returns:
10        ----------
11        labels: predicted cluster based on 
12        highest responsibility gamma.
13        
14        """
15        labels = np.zeros((X.shape[0], self.C))
16        
17        for c in range(self.C):
18            labels [:,c] = self.pi[c] * mvn.pdf(X, self.mu[c,:], self.sigma[c])
19        labels  = labels.argmax(1)
20        return labels 
21    
22    def predict_proba(self, X):
23        """Returns predicted labels
24        
25        Parameters:
26        -------------
27        X: N*d numpy array
28        
29        Returns:
30        ----------
31        labels: predicted cluster based on 
32        highest responsibility gamma.
33        
34        """
35        post_proba = np.zeros((X.shape[0], self.C))
36        
37        for c in range(self.C):
38            # Posterior Distribution using Bayes Rule
39            post_proba[:,c] = self.pi[c] * mvn.pdf(X, self.mu[c,:], self.sigma[c])
40    
41        return post_proba

1model = GMM(3, n_runs = 30)
2
3fitted_values = model.fit(Y_sklearn)
4predicted_values = model.predict(Y_sklearn)
5
6# # compute centers as point of highest density of distribution
7centers = np.zeros((3,2))
8for i in range(model.C):
9    density = mvn(cov=model.sigma[i], mean=model.mu[i]).logpdf(Y_sklearn)
10    centers[i, :] = Y_sklearn[np.argmax(density)]
11    
12plt.figure(figsize = (10,8))
13plt.scatter(Y_sklearn[:, 0], Y_sklearn[:, 1],c=predicted_values ,s=50, cmap='viridis', zorder=1)
14
15plt.scatter(centers[:, 0], centers[:, 1],c='black', s=300, alpha=0.5, zorder=2);
16
17w_factor = 0.2 / model.pi.max()
18for pos, covar, w in zip(model.mu, model.sigma, model.pi):
19    draw_ellipse(pos, covar, alpha = w)

1from sklearn.mixture import GaussianMixture
2sklearn_pca = PCA(n_components = 2)
3Y_sklearn = sklearn_pca.fit_transform(tf_idf_array)
4gmm = GaussianMixture(n_components=3, covariance_type='full').fit(Y_sklearn)
5prediction_gmm = gmm.predict(Y_sklearn)
6probs = gmm.predict_proba(Y_sklearn)
7
8centers = np.zeros((3,2))
9for i in range(3):
10    density = mvn(cov=gmm.covariances_[i], mean=gmm.means_[i]).logpdf(Y_sklearn)
11    centers[i, :] = Y_sklearn[np.argmax(density)]
12
13plt.figure(figsize = (10,8))
14plt.scatter(Y_sklearn[:, 0], Y_sklearn[:, 1],c=prediction_gmm ,s=50, cmap='viridis')
15plt.scatter(centers[:, 0], centers[:, 1],c='black', s=300, alpha=0.6);

1class GMM:
2    """ Gaussian Mixture Model
3    
4    Parameters
5    -----------
6        k: int , number of gaussian distributions
7        
8        seed: int, will be randomly set if None
9        
10        max_iter: int, number of iterations to run algorithm, default: 200
11        
12    Attributes
13    -----------
14       centroids: array, k, number_features
15       
16       cluster_labels: label for each data point
17       
18    """
19    def __init__(self, C, n_runs):
20        self.C = C # number of Guassians/clusters
21        self.n_runs = n_runs
22        
23    
24    def get_params(self):
25        return (self.mu, self.pi, self.sigma)
26    
27    
28        
29    def calculate_mean_covariance(self, X, prediction):
30        """Calculate means and covariance of different
31            clusters from k-means prediction
32        
33        Parameters:
34        ------------
35        prediction: cluster labels from k-means
36        
37        X: N*d numpy array data points 
38        
39        Returns:
40        -------------
41        intial_means: for E-step of EM algorithm
42        
43        intial_cov: for E-step of EM algorithm
44        
45        """
46        d = X.shape[1]
47        labels = np.unique(prediction)
48        self.initial_means = np.zeros((self.C, d))
49        self.initial_cov = np.zeros((self.C, d, d))
50        self.initial_pi = np.zeros(self.C)
51        
52        counter=0
53        for label in labels:
54            ids = np.where(prediction == label) # returns indices
55            self.initial_pi[counter] = len(ids[0]) / X.shape[0]
56            self.initial_means[counter,:] = np.mean(X[ids], axis = 0)
57            de_meaned = X[ids] - self.initial_means[counter,:]
58            Nk = X[ids].shape[0] # number of data points in current gaussian
59            self.initial_cov[counter,:, :] = np.dot(self.initial_pi[counter] * de_meaned.T, de_meaned) / Nk
60            counter+=1
61        assert np.sum(self.initial_pi) == 1    
62            
63        return (self.initial_means, self.initial_cov, self.initial_pi)
64    
65    
66    
67    def _initialise_parameters(self, X):
68        """Implement k-means to find starting
69            parameter values.
70            https://datascience.stackexchange.com/questions/11487/how-do-i-obtain-the-weight-and-variance-of-a-k-means-cluster
71
72        Parameters:
73        ------------
74        X: numpy array of data points
75        
76        Returns:
77        ----------
78        tuple containing initial means and covariance
79        
80        _initial_means: numpy array: (C*d)
81        
82        _initial_cov: numpy array: (C,d*d)
83        
84        
85        """
86        n_clusters = self.C
87        kmeans = KMeans(n_clusters= n_clusters, init="k-means++", max_iter=500, algorithm = 'auto')
88        fitted = kmeans.fit(X)
89        prediction = kmeans.predict(X)
90        self._initial_means, self._initial_cov, self._initial_pi = self.calculate_mean_covariance(X, prediction)
91        
92        
93        return (self._initial_means, self._initial_cov, self._initial_pi)
94            
95        
96        
97    def _e_step(self, X, pi, mu, sigma):
98        """Performs E-step on GMM model
99
100        Parameters:
101        ------------
102        X: (N x d), data points, m: no of features
103        pi: (C), weights of mixture components
104        mu: (C x d), mixture component means
105        sigma: (C x d x d), mixture component covariance matrices
106
107        Returns:
108        ----------
109        gamma: (N x C), probabilities of clusters for objects
110        """
111        N = X.shape[0] 
112        self.gamma = np.zeros((N, self.C))
113
114        const_c = np.zeros(self.C)
115        
116        
117        self.mu = self.mu if self._initial_means is None else self._initial_means
118        self.pi = self.pi if self._initial_pi is None else self._initial_pi
119        self.sigma = self.sigma if self._initial_cov is None else self._initial_cov
120
121        for c in range(self.C):
122            # Posterior Distribution using Bayes Rule
123            self.gamma[:,c] = self.pi[c] * mvn.pdf(X, self.mu[c,:], self.sigma[c])
124
125        # normalize across columns to make a valid probability
126        gamma_norm = np.sum(self.gamma, axis=1)[:,np.newaxis]
127        self.gamma /= gamma_norm
128
129        return self.gamma
130    
131    
132    def _m_step(self, X, gamma):
133        """Performs M-step of the GMM
134        We need to update our priors, our means
135        and our covariance matrix.
136
137        Parameters:
138        -----------
139        X: (N x d), data 
140        gamma: (N x C), posterior distribution of lower bound 
141
142        Returns:
143        ---------
144        pi: (C)
145        mu: (C x d)
146        sigma: (C x d x d)
147        """
148        N = X.shape[0] # number of objects
149        C = self.gamma.shape[1] # number of clusters
150        d = X.shape[1] # dimension of each object
151
152        # responsibilities for each gaussian
153        self.pi = np.mean(self.gamma, axis = 0)
154
155        self.mu = np.dot(self.gamma.T, X) / np.sum(self.gamma, axis = 0)[:,np.newaxis]
156
157        for c in range(C):
158            x = X - self.mu[c, :] # (N x d)
159            
160            gamma_diag = np.diag(self.gamma[:,c])
161            x_mu = np.matrix(x)
162            gamma_diag = np.matrix(gamma_diag)
163
164            sigma_c = x.T * gamma_diag * x
165            self.sigma[c,:,:]=(sigma_c) / np.sum(self.gamma, axis = 0)[:,np.newaxis][c]
166
167        return self.pi, self.mu, self.sigma
168    
169    
170    def _compute_loss_function(self, X, pi, mu, sigma):
171        """Computes lower bound loss function
172        
173        Parameters:
174        -----------
175        X: (N x d), data 
176        
177        Returns:
178        ---------
179        pi: (C)
180        mu: (C x d)
181        sigma: (C x d x d)
182        """
183        N = X.shape[0]
184        C = self.gamma.shape[1]
185        self.loss = np.zeros((N, C))
186
187        for c in range(C):
188            dist = mvn(self.mu[c], self.sigma[c],allow_singular=True)
189            self.loss[:,c] = self.gamma[:,c] * (np.log(self.pi[c]+0.00001)+dist.logpdf(X)-np.log(self.gamma[:,c]+0.000001))
190        self.loss = np.sum(self.loss)
191        return self.loss
192    
193    
194    
195    def fit(self, X):
196        """Compute the E-step and M-step and
197            Calculates the lowerbound
198        
199        Parameters:
200        -----------
201        X: (N x d), data 
202        
203        Returns:
204        ----------
205        instance of GMM
206        
207        """
208        
209        d = X.shape[1]
210        self.mu, self.sigma, self.pi =  self._initialise_parameters(X)
211        
212        try:
213            for run in range(self.n_runs):  
214                self.gamma  = self._e_step(X, self.mu, self.pi, self.sigma)
215                self.pi, self.mu, self.sigma = self._m_step(X, self.gamma)
216                loss = self._compute_loss_function(X, self.pi, self.mu, self.sigma)
217                
218                if run % 10 == 0:
219                    print("Iteration: %d Loss: %0.6f" %(run, loss))
220
221        
222        except Exception as e:
223            print(e)
224            
225        
226        return self
227    
228    
229    
230    
231    def predict(self, X):
232        """Returns predicted labels using Bayes Rule to
233        Calculate the posterior distribution
234        
235        Parameters:
236        -------------
237        X: ?*d numpy array
238        
239        Returns:
240        ----------
241        labels: predicted cluster based on 
242        highest responsibility gamma.
243        
244        """
245        labels = np.zeros((X.shape[0], self.C))
246        
247        for c in range(self.C):
248            labels [:,c] = self.pi[c] * mvn.pdf(X, self.mu[c,:], self.sigma[c])
249        labels  = labels .argmax(1)
250        return labels 
251    
252    def predict_proba(self, X):
253        """Returns predicted labels
254        
255        Parameters:
256        -------------
257        X: N*d numpy array
258        
259        Returns:
260        ----------
261        labels: predicted cluster based on 
262        highest responsibility gamma.
263        
264        """
265        post_proba = np.zeros((X.shape[0], self.C))
266        
267        for c in range(self.C):
268            # Posterior Distribution using Bayes Rule, try and vectorise
269            post_proba[:,c] = self.pi[c] * mvn.pdf(X, self.mu[c,:], self.sigma[c])
270    
271        return post_proba