# How does scikit-learn calculate a data point's probability of belong to normal distribution?

In GMM calculating, we need to calculate the probability of data point $X$ belong to the $kth$ Gaussian Distribution $\mathcal{N}(X_n|\mu_k,\Sigma_k)$.

I have read How to calculate the probability of a data point belonging to a multivariate normal distribution? but still having trouble understanding how scikit-learn calculated it.

The value returned is the log probability of $\mathcal{N}(X_n|\mu_k,\Sigma_k)$. Would someone please explain what does this return value mean?

The following code is from sklearn.mixture.GaussianMixture:

def _estimate_log_gaussian_prob(X, means, precisions_chol, covariance_type):
"""Estimate the log Gaussian probability.

Parameters
----------
X : array-like, shape (n_samples, n_features)

means : array-like, shape (n_components, n_features)

precisions_chol : array-like,
Cholesky decompositions of the precision matrices.
'full' : shape of (n_components, n_features, n_features)
'tied' : shape of (n_features, n_features)
'diag' : shape of (n_components, n_features)
'spherical' : shape of (n_components,)

covariance_type : {'full', 'tied', 'diag', 'spherical'}

Returns
-------
log_prob : array, shape (n_samples, n_components)
"""
n_samples, n_features = X.shape
n_components, _ = means.shape
# det(precision_chol) is half of det(precision)
log_det = _compute_log_det_cholesky(
precisions_chol, covariance_type, n_features)

if covariance_type == 'full':
log_prob = np.empty((n_samples, n_components))
for k, (mu, prec_chol) in enumerate(zip(means, precisions_chol)):
y = np.dot(X, prec_chol) - np.dot(mu, prec_chol)
log_prob[:, k] = np.sum(np.square(y), axis=1)

elif covariance_type == 'tied':
log_prob = np.empty((n_samples, n_components))
for k, mu in enumerate(means):
y = np.dot(X, precisions_chol) - np.dot(mu, precisions_chol)
log_prob[:, k] = np.sum(np.square(y), axis=1)

elif covariance_type == 'diag':
precisions = precisions_chol ** 2
log_prob = (np.sum((means ** 2 * precisions), 1) -
2. * np.dot(X, (means * precisions).T) +
np.dot(X ** 2, precisions.T))

elif covariance_type == 'spherical':
precisions = precisions_chol ** 2
log_prob = (np.sum(means ** 2, 1) * precisions -
2 * np.dot(X, means.T * precisions) +
np.outer(row_norms(X, squared=True), precisions))
return -.5 * (n_features * np.log(2 * np.pi) + log_prob) + log_det


As you have stated in the question the log probability of the multivariate gaussian is as follows. I'll take a single component to simplify notation. $$\log{\mathcal{N}(x | \mu, \Sigma)} = -\frac{1}{2}K\log{2\pi} - \frac{1}{2}\log{|\Sigma|} - \frac{1}{2} (x - \mu)^T\Sigma^{-1}(x-\mu) \tag{1}$$

The Scipy algorithm is using the Cholesky decomposition of the precision matrices. The precision is the inverse of the covariance, i.e., $$\Lambda = \Sigma^{-1}$$ Let $$L$$ denote the Cholesky decomposition of a particular precision matrix, so that, $$\Lambda= LL^T$$

Computing the log_prob

Starting from the Mahalanobis term in Equation (1), \begin{align} ( x - \mu )^T \Sigma^{-1} (x - \mu ) =& ( x - \mu )^T \Lambda (x - \mu ) \\ =& ( x - \mu )^T LL^T (x - \mu ) \\ =& \|L^T (x - \mu )\|^2 \end{align}

The algorithm first computes $$y = L^T (x - \mu_i )$$ using the dot function, and then takes the norm squared, i.e. $$log\_prob = \|y\|^2$$. It then needs to multiply by a factor of $$-\frac{1}{2}$$, which happens in the return statement.

Computing the log_det

Starting from the log determinant term in equation (1), \begin{align} -\frac{1}{2}\log{|\Sigma|} =& -\frac{1}{2} \log{|\Lambda^{-1}|} \\ =& \frac{1}{2}\log{|\Lambda|} \\ =& \frac{1}{2}\log{|LL^T|} \\ =& \log{|L|} \end{align}

In the return statement, $$log\_det = \det{L}$$, and so there is no need to multiply by $$-\frac{1}{2}$$.

The whole computation

In summary, the final computation in Scipy is, $$\log{\mathcal{N}(x | \mu, \Sigma)} = -\frac{1}{2} \left( K\log{2\pi} + \|L^T(x - \mu)\|^2 \right) +\log{|L|}$$