slow_computation_gaussian_mixture.py

# Compute the probability of each component in the mixture
def _computeProbabilities(self, limits):
    probabilities = np.zeros(numComponents) # numComponents refers to the number of components in this mixture
    for component in range(numComponents):
		# gmm is a GaussianMixture from sklearn.mixture
        mean = gmm.means_[component]
        cov = gmm.covariances_[component]
        weight = gmm.weights_[component]
        probabilities[component] = constrainedGaussianDensity(mean, cov, limits) * weight
    return probabilities

# Calculates total density within a constrained rectangle of a Gaussian with 
# specified mean and (diagonal) covariance matrix
#
#
# params/returns:
#  mean: np.array([nCols]), (multivariate) mean of the distribution
#  cov: np.array([nCols]), (multivariate) diagonal elements of the covariance matrix of the distribution, assumed diagonal
#  limits: [[float | np.inf | -np.inf]], a list of lists, each of which specifies an upper and lower bound
#    of the constrained region for that dimension
#  return: float
def constrainedGaussianDensity(mean, cov, limits):
    nCols = np.size(mean)
    totalDensity = 1
    for d in range(nCols):
        limit = limits[d]
        std = np.sqrt(cov[d])
        totalDensity = totalDensity * integrateTruncated(mean[d], std, limit[0], limit[1])
    return totalDensity

# Calculates integral under a 1-dimensional Gaussian with specified mean, 
# standard deviation, and between the two ranges specified
# parameters/returns:
#  mean : float, mean of distribution
#  std: float, standard deviation of distribution
#  lower : float | -np.inf, lower limit of integral 
#  upper : float | np.inf, upper limit of integral 
#  return : float, area between the two limits
def integrateTruncated(mean, std, lower, upper):
    rv = norm(loc=mean, scale=std)
    if rv.cdf(upper) - rv.cdf(lower) == 0:
        log(("equals zero", upper, lower, mean, std), DEBUG)
    return rv.cdf(upper) - rv.cdf(lower)