'Compute data covariance with the FactorAnalysis model.
``cov = components_.T * components_ + diag(noise_variance)``
Returns
cov : array, shape (n_features, n_features)
Estimated covariance of data.'

def get_covariance(self):

    check_is_fitted(self, 'components_')
    cov = np.dot(self.components_.T, self.components_)
    cov.flat[::(len(cov) + 1)] += self.noise_variance_
    return cov
