Python extmath.pinvh函数代码示例

 def fit(self, balance_param=0.5, sparsity_param=0.01, verbose=False):
   '''
   balance_param: trades off between sparsity and M0 prior
   sparsity_param: trades off between optimizer and sparseness (see graph_lasso)
   '''
   P = pinvh(self.M) + balance_param * self.loss_matrix
   emp_cov = pinvh(P)
   # hack: ensure positive semidefinite
   emp_cov = emp_cov.T.dot(emp_cov)
   self.M, _ = graph_lasso(emp_cov, sparsity_param, verbose=verbose)

 def fit(self, X, W, verbose=False):
   """
   X: data matrix, (n x d)
   W: connectivity graph, (n x n). +1 for positive pairs, -1 for negative.
   """
   self._prepare_inputs(X, W)
   P = pinvh(self.M) + self.balance_param * self.loss_matrix
   emp_cov = pinvh(P)
   # hack: ensure positive semidefinite
   emp_cov = emp_cov.T.dot(emp_cov)
   self.M, _ = graph_lasso(emp_cov, self.sparsity_param, verbose=verbose)
   return self

 def fit(self, X, W):
   """
   X: data matrix, (n x d)
       each row corresponds to a single instance
   W: connectivity graph, (n x n). +1 for positive pairs, -1 for negative.
   """
   self._prepare_inputs(X, W)
   P = pinvh(self.M) + self.params['balance_param'] * self.loss_matrix
   emp_cov = pinvh(P)
   # hack: ensure positive semidefinite
   emp_cov = emp_cov.T.dot(emp_cov)
   self.M, _ = graph_lasso(emp_cov, self.params['sparsity_param'],
                           verbose=self.params['verbose'])
   return self

    def correct_covariance(self, data, method=None):
        """Apply a correction to raw Minimum Covariance Determinant estimates.

        Correction using the empirical correction factor suggested
        by Rousseeuw and Van Driessen in [Rouseeuw1984]_.

        Parameters
        ----------
        data: array-like, shape (n_samples, n_features)
          The data matrix, with p features and n samples.
          The data set must be the one which was used to compute
          the raw estimates.

        Returns
        -------
        covariance_corrected: array-like, shape (n_features, n_features)
          Corrected robust covariance estimate.

        """
        if method is "empirical":
            X_c = data - self.raw_location_
            dist = np.sum(
                np.dot(X_c, pinvh(self.raw_covariance_)) * X_c, 1)
            correction = np.median(dist) / sp.stats.chi2(
                data.shape[1]).isf(0.5)
            covariance_corrected = self.raw_covariance_ * correction
        elif method is "theoretical":
            n, p = data.shape
            c = sp.stats.chi2(p + 2).cdf(sp.stats.chi2(p).ppf(self.h)) / self.h
            covariance_corrected = self.raw_covariance_ * c
        else:
            covariance_corrected = self.raw_covariance_
        self._set_covariance(covariance_corrected)

        return covariance_corrected

 def _posterior_dist(self,A,beta,XX,XY,full_covar=False):
     '''
     Calculates mean and covariance matrix of posterior distribution
     of coefficients.
     '''
     # compute precision matrix for active features
     Sinv = beta * XX
     np.fill_diagonal(Sinv, np.diag(Sinv) + A)
     cholesky = True
     # try cholesky, if it fails go back to pinvh
     try:
         # find posterior mean : R*R.T*mean = beta*X.T*Y
         # solve(R*z = beta*X.T*Y) => find z => solve(R.T*mean = z) => find mean
         R    = np.linalg.cholesky(Sinv)
         Z    = solve_triangular(R,beta*XY, check_finite=False, lower = True)
         Mn   = solve_triangular(R.T,Z, check_finite=False, lower = False)
         
         # invert lower triangular matrix from cholesky decomposition
         Ri   = solve_triangular(R,np.eye(A.shape[0]), check_finite=False, lower=True)
         if full_covar:
             Sn   = np.dot(Ri.T,Ri)
             return Mn,Sn,cholesky
         else:
             return Mn,Ri,cholesky
     except LinAlgError:
         cholesky = False
         Sn   = pinvh(Sinv)
         Mn   = beta*np.dot(Sinv,XY)
         return Mn, Sn, cholesky

def mutual_incoherence(X_relevant, X_irelevant):
    """Mutual incoherence, as defined by formula (26a) of [Wainwright2006].
    """
    projector = np.dot(
                    np.dot(X_irelevant.T, X_relevant),
                    pinvh(np.dot(X_relevant.T, X_relevant))
                    )
    return np.max(np.abs(projector).sum(axis=1))

def _log_multivariate_normal_density_tied(X, means, covars):
    """Compute Gaussian log-density at X for a tied model"""
    n_samples, n_dim = X.shape
    icv = pinvh(covars)
    lpr = -0.5 * (n_dim * np.log(2 * np.pi) + np.log(linalg.det(covars) + 0.1)
                  + np.sum(X * np.dot(X, icv), 1)[:, np.newaxis]
                  - 2 * np.dot(np.dot(X, icv), means.T)
                  + np.sum(means * np.dot(means, icv), 1))
    return lpr

def test_pinvh_nonpositive():
    a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=np.float64)
    a = np.dot(a, a.T)
    u, s, vt = np.linalg.svd(a)
    s[0] *= -1
    a = np.dot(u * s, vt)  # a is now symmetric non-positive and singular
    a_pinv = pinv2(a)
    a_pinvh = pinvh(a)
    assert_almost_equal(a_pinv, a_pinvh)

 def objective_function(self, data, location, covariance):
     """Objective function minimized at each step of the MCD algorithm.
     """
     precision = pinvh(covariance)
     det = fast_logdet(precision)
     trace = np.trace(
         np.dot(empirical_covariance(data - location, assume_centered=True),
                precision))
     pen = self.shrinkage * np.trace(precision)
     return -det + trace + pen

 def _prepare_inputs(self, X, W):
   self.X_ = X = check_array(X)
   W = check_array(W, accept_sparse=True)
   # set up prior M
   if self.use_cov:
     self.M_ = pinvh(np.cov(X, rowvar = False))
   else:
     self.M_ = np.identity(X.shape[1])
   L = laplacian(W, normed=False)
   return X.T.dot(L.dot(X))

    def _condition(self, i1, i2, X):
        cov_12 = self.covariance[np.ix_(i1, i2)]
        cov_11 = self.covariance[np.ix_(i1, i1)]
        cov_22 = self.covariance[np.ix_(i2, i2)]
        prec_22 = pinvh(cov_22)
        regression_coeffs = cov_12.dot(prec_22)

        if X.ndim == 2:
            mean = self.mean[i1] + regression_coeffs.dot(
                (X - self.mean[i2]).T).T
        elif X.ndim == 1:
            mean = self.mean[i1] + regression_coeffs.dot(X - self.mean[i2])
        else:
            raise ValueError("%d dimensions are not allowed for X!" % X.ndim)
        covariance = cov_11 - regression_coeffs.dot(cov_12.T)
        return mean, covariance

    def fit(self, X, y=None):
        """Fits a Minimum Covariance Determinant with the FastMCD algorithm.

        Parameters
        ----------
        X: array-like, shape = [n_samples, n_features]
          Training data, where n_samples is the number of samples
          and n_features is the number of features.
        y: not used, present for API consistence purpose.

        Returns
        -------
        self: object
          Returns self.

        """
        n_samples, n_features = X.shape
        # check that the empirical covariance is full rank
        if (linalg.svdvals(np.dot(X.T, X)) > 1e-8).sum() != n_features:
            warnings.warn("The covariance matrix associated to your dataset "
                          "is not full rank")
        # compute and store raw estimates
        raw_location, raw_covariance, raw_support = fast_mcd(
            X, objective_function=self.objective_function,
            h=self.h, cov_computation_method=self._nonrobust_covariance)
        if self.h is None:
            self.h = int(np.ceil(0.5 * (n_samples + n_features + 1))) \
                / float(n_samples)
        if self.assume_centered:
            raw_location = np.zeros(n_features)
            raw_covariance = self._nonrobust_covariance(
                X[raw_support], assume_centered=True)
        # get precision matrix in an optimized way
        precision = pinvh(raw_covariance)
        raw_dist = np.sum(np.dot(X, precision) * X, 1)
        self.raw_location_ = raw_location
        self.raw_covariance_ = raw_covariance
        self.raw_support_ = raw_support
        self.location_ = raw_location
        self.support_ = raw_support
        self.dist_ = raw_dist
        # obtain consistency at normal models
        self.correct_covariance(X)

        return self

def launch_rmcdl1_on_dataset(n_samples, n_features, n_outliers):

    rand_gen = np.random.RandomState(0)
    data = rand_gen.randn(n_samples, n_features)
    # add some outliers
    outliers_index = rand_gen.permutation(n_samples)[:n_outliers]
    outliers_offset = 10. * \
        (rand_gen.randint(2, size=(n_outliers, n_features)) - 0.5)
    data[outliers_index] += outliers_offset
    inliers_mask = np.ones(n_samples).astype(bool)
    inliers_mask[outliers_index] = False

    # compute RMCD by fitting an object
    rmcd_fit = RMCDl1().fit(data)
    T = rmcd_fit.location_
    S = rmcd_fit.covariance_
    # compare with the true location and precision
    error_location = np.mean(T ** 2)
    assert(error_location < 1.)
    error_cov = np.mean((np.eye(n_features) - pinvh(S)) ** 2)
    assert(error_cov < 1.)

 def _posterior_dist(self,X,y,A,intercept_prior):
     '''
     Uses Laplace approximation for calculating posterior distribution
     '''
     if self.solver == 'lbfgs_b':
         f  = lambda w: _logistic_cost_grad(X,y,w,A,intercept_prior)
         w_init  = np.random.random(X.shape[1])
         Mn      = fmin_l_bfgs_b(f, x0 = w_init, pgtol = self.tol_solver,
                                 maxiter = self.n_iter_solver)[0]
         Xm      = np.dot(X,Mn)
         s       = expit(Xm)
         B       = logistic._pdf(Xm) # avoids underflow
         S       = np.dot(X.T*B,X)
         np.fill_diagonal(S, np.diag(S) + A)
         t_hat   = Xm + (y - s) / B
         Sn      = pinvh(S)
     elif self.solver == 'newton_cg':
         # TODO: Implement Newton-CG
         raise NotImplementedError(('Newton Conjugate Gradient optimizer '
                                    'is not currently supported'))
     return [Mn,Sn,B,t_hat]

  def fit(self, X, W):
    """Learn the SDML model.

    Parameters
    ----------
    X : array-like, shape (n, d)
        data matrix, where each row corresponds to a single instance
    W : array-like, shape (n, n)
        connectivity graph, with +1 for positive pairs and -1 for negative

    Returns
    -------
    self : object
        Returns the instance.
    """
    loss_matrix = self._prepare_inputs(X, W)
    P = self.M_ + self.balance_param * loss_matrix
    emp_cov = pinvh(P)
    # hack: ensure positive semidefinite
    emp_cov = emp_cov.T.dot(emp_cov)
    _, self.M_ = graph_lasso(emp_cov, self.sparsity_param, verbose=self.verbose)
    return self