'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]
    dg = (digamma(self.gamma_) - digamma(np.sum(self.gamma_)))
    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 + dg)
    z = log_normalize(z, axis=(-1))
    bound = np.sum((z * p), axis=(-1))
    return (bound, z)
