'Estimate model parameters with the variational
algorithm.
For a full derivation and description of the algorithm see
doc/modules/dp-derivation.rst
or
http://scikit-learn.org/stable/modules/dp-derivation.html
A initialization step is performed before entering the em
algorithm. If you want to avoid this step, set the keyword
argument init_params to the empty string \'\' when creating
the object. Likewise, if you would like just to do an
initialization, set n_iter=0.
Parameters
X : array_like, shape (n, n_features)
List of n_features-dimensional data points.  Each row
corresponds to a single data point.
Returns
responsibilities : array, shape (n_samples, n_components)
Posterior probabilities of each mixture component for each
observation.'

def _fit(self, X, y=None):

    self.random_state_ = check_random_state(self.random_state)
    X = check_array(X)
    if (X.ndim == 1):
        X = X[:, np.newaxis]
    (n_samples, n_features) = X.shape
    z = np.ones((n_samples, self.n_components))
    z /= self.n_components
    self._initial_bound = (((-0.5) * n_features) * np.log((2 * np.pi)))
    self._initial_bound -= np.log(((2 * np.pi) * np.e))
    if ((self.init_params != '') or (not hasattr(self, 'gamma_'))):
        self._initialize_gamma()
    if (('m' in self.init_params) or (not hasattr(self, 'means_'))):
        self.means_ = cluster.KMeans(n_clusters=self.n_components, random_state=self.random_state_).fit(X).cluster_centers_[::(-1)]
    if (('w' in self.init_params) or (not hasattr(self, 'weights_'))):
        self.weights_ = np.tile((1.0 / self.n_components), self.n_components)
    if (('c' in self.init_params) or (not hasattr(self, 'precs_'))):
        if (self.covariance_type == 'spherical'):
            self.dof_ = np.ones(self.n_components)
            self.scale_ = np.ones(self.n_components)
            self.precs_ = np.ones((self.n_components, n_features))
            self.bound_prec_ = ((0.5 * n_features) * (digamma(self.dof_) - np.log(self.scale_)))
        elif (self.covariance_type == 'diag'):
            self.dof_ = (1 + (0.5 * n_features))
            self.dof_ *= np.ones((self.n_components, n_features))
            self.scale_ = np.ones((self.n_components, n_features))
            self.precs_ = np.ones((self.n_components, n_features))
            self.bound_prec_ = (0.5 * np.sum((digamma(self.dof_) - np.log(self.scale_)), 1))
            self.bound_prec_ -= (0.5 * np.sum(self.precs_, 1))
        elif (self.covariance_type == 'tied'):
            self.dof_ = 1.0
            self.scale_ = np.identity(n_features)
            self.precs_ = np.identity(n_features)
            self.det_scale_ = 1.0
            self.bound_prec_ = (0.5 * wishart_log_det(self.dof_, self.scale_, self.det_scale_, n_features))
            self.bound_prec_ -= ((0.5 * self.dof_) * np.trace(self.scale_))
        elif (self.covariance_type == 'full'):
            self.dof_ = ((1 + self.n_components) + n_samples)
            self.dof_ *= np.ones(self.n_components)
            self.scale_ = [(2 * np.identity(n_features)) for _ in range(self.n_components)]
            self.precs_ = [np.identity(n_features) for _ in range(self.n_components)]
            self.det_scale_ = np.ones(self.n_components)
            self.bound_prec_ = np.zeros(self.n_components)
            for k in range(self.n_components):
                self.bound_prec_[k] = wishart_log_det(self.dof_[k], self.scale_[k], self.det_scale_[k], n_features)
                self.bound_prec_[k] -= (self.dof_[k] * np.trace(self.scale_[k]))
            self.bound_prec_ *= 0.5
    current_log_likelihood = None
    self.converged_ = False
    for i in range(self.n_iter):
        prev_log_likelihood = current_log_likelihood
        (curr_logprob, z) = self.score_samples(X)
        current_log_likelihood = (curr_logprob.mean() + (self._logprior(z) / n_samples))
        if (prev_log_likelihood is not None):
            change = abs((current_log_likelihood - prev_log_likelihood))
            if (change < self.tol):
                self.converged_ = True
                break
        self._do_mstep(X, z, self.params)
    if (self.n_iter == 0):
        z = np.zeros((X.shape[0], self.n_components))
    self._set_weights()
    return z
