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

 A: Please check this link: 
https://math.stackexchange.com/a/2673224/707024
It's really helpful
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|}
$$
