'Return the likelihood of the data under the model.
Compute the bound on log probability of X under the model
and return the posterior distribution (responsibilities) of
each mixture component for each element of X.
This is done by computing the parameters for the mean-field of
z for each observation.
Parameters
X : array_like, shape (n_samples, n_features)
List of n_features-dimensional data points.  Each row
corresponds to a single data point.
Returns
logprob : array_like, shape (n_samples,)
Log probabilities of each data point in X
responsibilities : array_like, shape (n_samples, n_components)
Posterior probabilities of each mixture component for each
observation'

def score_samples(self, X):

    check_is_fitted(self, 'gamma_')
    X = check_array(X)
    if (X.ndim == 1):
        X = X[:, np.newaxis]
    z = np.zeros((X.shape[0], self.n_components))
    sd = digamma((self.gamma_.T[1] + self.gamma_.T[2]))
    dgamma1 = (digamma(self.gamma_.T[1]) - sd)
    dgamma2 = np.zeros(self.n_components)
    dgamma2[0] = (digamma(self.gamma_[(0, 2)]) - digamma((self.gamma_[(0, 1)] + self.gamma_[(0, 2)])))
    for j in range(1, self.n_components):
        dgamma2[j] = (dgamma2[(j - 1)] + digamma(self.gamma_[((j - 1), 2)]))
        dgamma2[j] -= sd[(j - 1)]
    dgamma = (dgamma1 + dgamma2)
    del dgamma1, dgamma2, sd
    if (self.covariance_type not in ['full', 'tied', 'diag', 'spherical']):
        raise NotImplementedError(('This    ctype    is    not    implemented:    %s' % self.covariance_type))
    p = _bound_state_log_lik(X, (self._initial_bound + self.bound_prec_), self.precs_, self.means_, self.covariance_type)
    z = (p + dgamma)
    z = log_normalize(z, axis=(-1))
    bound = np.sum((z * p), axis=(-1))
    return (bound, z)
