Getting the PDF from the Gausian Mixture Model in sklearn

Viewed 376

I have fit a Gaussian Mixture Model (GMM) to a data series that I have. Using GMM, I am trying to get the probabilities of another vector, element-wise. Matlab achieves this with the following lines of code.

a = reshape(0:1:15, 14, 1);
gm = fitgmdist(a, 13);  % 13 specifies the number of components (means and covs for example) in the fit model

% Testing with new data
b = reshape(-5:1:5, 11, 1);
pdf(gm, b) 

    0.0000
    0.0000
    0.0000
    0.0000
    0.0018
    0.0643
    0.0658
    0.0671
    0.0666
    0.0662
    0.0672

This is expected because the negative values -5 through 0 are not present in the data provided while fitting and therefore provide near zero values.

I am trying to replicate this using python with sklearn. Following is what I have achieved so far.

from sklearn.mixture import GaussianMixture
import numpy as np

gm = GaussianMixture(n_components=13).fit(np.arange(16).reshape(-1, 1))

# Generate test data
b = np.arange(-5, 6)[:, None]
prob = gm.predict_proba(b).tolist()

"""
prob=
[[8.539939840944505e-152, 0.0, 1.9638856033086253e-68, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[2.238593143299414e-141, 0.0, 3.166463050557315e-63, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[5.868073258732079e-131, 0.0, 5.106947259415683e-58, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[1.5382109014584666e-120, 0.0, 8.239047963148164e-53, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[4.0321459413005606e-110, 0.0, 1.3296012030592655e-47, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[8.790691487706948e-139, 0.0, 1.7850932827696813e-81, 8.316994953954272e-40, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[6.040819196361928e-118, 0.0, 7.556403579572576e-66, 2.180313173877784e-29, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0], 
[7.616317001277741e-99, 0.0, 5.870491452246584e-52, 1.0486913731065672e-20, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[1.761856615767359e-81, 0.9999999999999076, 8.370258307836422e-40, 9.254498458447064e-14, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 
[7.477776611058069e-66, 0.0, 2.190323924113167e-29, 1.4984230843298664e-08, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9999999850157694, 0.0, 0.0], 
[5.823053603579932e-52, 0.0, 1.0519204995101235e-20, 4.4513567127132954e-05, 0.0, 0.0, 0.9999554864328727, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
"""

The array prob returns a matrix in which each row signifies the probability of belonging to one of the n_components number of classes (implying sum along row = 1). This is however not the expected result. I would want to know the probability of each element of b being generated from the fit gm model, like Matlab does.

How can I achieve this through Python? Thanks.

1 Answers

predict_proba finds the probability that the individual data row is assigned to each of component the GMM is based on (in your case, 13 components). For more, see Jake VanderPlas' excellent chapter.

To get the pdf, you'll need to use the score_samples function-which returns the weighted log probability that each row would be drawn from the fitted GMM:

b = np.arange(-5, 6)[:, None]
log_probs = gm.score_samples(b)
sum_of_scores = (np.exp(log_probs)).sum()
probs = np.exp(log_probs) / sum_of_scores
print('Probs sum: ', probs.sum()) # confirm sums to 1
print('pdf: ', probs.tolist())

Note that in this toy example, due to the non-deterministic nature of the K-Means clustering at the initial part of the GMM fit algorithm, you most likely will not reproduce the exact results from matlab, but if you do have a more realistic, diverse, well sampled dataset to cross-check against the MATLAB implementation, they will agree much better.

Related