Python linalg.lstsq方法代码示例

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def infer(self, Ms, ys, scale_factors=None):
        ''' Either:
            1) Ms is a single M and ys is a single y 
               (scale_factors ignored) or
            2) Ms and ys are lists of M matrices and y vectors
               and scale_factors is a list of the same length.
        '''
        A, y = _apply_scales(Ms, ys, scale_factors)

        if self.method == 'standard':
            assert self.l2_reg == 0, 'l2 reg not supported with method=standard'
            (x_est, _, rank, _) = linalg.lstsq(A.dense_matrix(), y, lapack_driver='gelsy')
        elif self.method == 'lsmr':
            res = lsmr(A, y, atol=0, btol=0, damp=self.l2_reg)
            x_est = res[0]
        elif self.method == 'lsqr':
            res = lsqr(A, y, atol=0, btol=0, damp=self.l2_reg)
            x_est = res[0]

        return x_est 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def mtx_updated_G(phi_recon, M, mtx_amp2visi_ri, mtx_fri2visi_ri):
    """
    Update the linear transformation matrix that links the FRI sequence to the
    visibilities by using the reconstructed Dirac locations.
    :param phi_recon: the reconstructed Dirac locations (azimuths)
    :param M: the Fourier series expansion is between -M to M
    :param p_mic_x: a vector that contains microphones' x-coordinates
    :param p_mic_y: a vector that contains microphones' y-coordinates
    :param mtx_freq2visi: the linear mapping from Fourier series to visibilities
    :return:
    """
    L = 2 * M + 1
    ms_half = np.reshape(np.arange(-M, 1, step=1), (-1, 1), order='F')
    phi_recon = np.reshape(phi_recon, (1, -1), order='F')
    mtx_amp2freq = np.exp(-1j * ms_half * phi_recon)  # size: (M + 1) x K
    mtx_amp2freq_ri = np.vstack((mtx_amp2freq.real, mtx_amp2freq.imag[:-1, :]))  # size: (2M + 1) x K
    mtx_fri2amp_ri = linalg.lstsq(mtx_amp2freq_ri, np.eye(L))[0]
    # projection mtx_freq2visi to the null space of mtx_fri2amp
    mtx_null_proj = np.eye(L) - np.dot(mtx_fri2amp_ri.T,
                                       linalg.lstsq(mtx_fri2amp_ri.T, np.eye(L))[0])
    G_updated = np.dot(mtx_amp2visi_ri, mtx_fri2amp_ri) + \
                np.dot(mtx_fri2visi_ri, mtx_null_proj)
    return G_updated 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def _update_means(self, X, z):
        """Update the variational distributions for the means"""
        n_features = X.shape[1]
        for k in range(self.n_components):
            if self.covariance_type in ['spherical', 'diag']:
                num = np.sum(z.T[k].reshape((-1, 1)) * X, axis=0)
                num *= self.precs_[k]
                den = 1. + self.precs_[k] * np.sum(z.T[k])
                self.means_[k] = num / den
            elif self.covariance_type in ['tied', 'full']:
                if self.covariance_type == 'tied':
                    cov = self.precs_
                else:
                    cov = self.precs_[k]
                den = np.identity(n_features) + cov * np.sum(z.T[k])
                num = np.sum(z.T[k].reshape((-1, 1)) * X, axis=0)
                num = np.dot(cov, num)
                self.means_[k] = linalg.lstsq(den, num)[0] 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def alt_orbital_correction(self, ifg, deg, offset):
        data = ifg.phase_data.reshape(ifg.num_cells)
        dm = get_design_matrix(ifg, deg, offset)[~isnan(data)]
        fd = data[~isnan(data)].reshape((dm.shape[0], 1))

        dmt = dm.T
        invNbb = inv(dmt.dot(dm))
        orbparams = invNbb.dot(dmt.dot(fd))
        alt_params = lstsq(dm, fd)[0]
        # FIXME: precision
        assert_array_almost_equal(orbparams, alt_params, decimal=2)

        dm2 = get_design_matrix(ifg, deg, offset)

        if offset:
            fullorb = np.reshape(np.dot(dm2[:, :-1], orbparams[:-1]),
                             ifg.phase_data.shape)
        else:
            fullorb = np.reshape(np.dot(dm2, orbparams), ifg.phase_data.shape)

        offset_removal = nanmedian(
            np.reshape(ifg.phase_data - fullorb, (1, -1)))
        fwd_correction = fullorb - offset_removal
        # ifg.phase_data -= (fullorb - offset_removal)
        return ifg.phase_data - fwd_correction 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def mtx_updated_G_multiband(phi_recon, M, mtx_amp2visi_ri,
                            mtx_fri2visi_ri, num_bands):
    """
    Update the linear transformation matrix that links the FRI sequence to the
    visibilities by using the reconstructed Dirac locations.
    :param phi_recon: the reconstructed Dirac locations (azimuths)
    :param M: the Fourier series expansion is between -M to M
    :param p_mic_x: a vector that contains microphones' x-coordinates
    :param p_mic_y: a vector that contains microphones' y-coordinates
    :param mtx_freq2visi: the linear mapping from Fourier series to visibilities
    :return:
    """
    L = 2 * M + 1
    ms_half = np.reshape(np.arange(-M, 1, step=1), (-1, 1), order='F')
    phi_recon = np.reshape(phi_recon, (1, -1), order='F')
    mtx_amp2freq = np.exp(-1j * ms_half * phi_recon)  # size: (M + 1) x K
    mtx_amp2freq_ri = np.vstack((mtx_amp2freq.real, mtx_amp2freq.imag[:-1, :]))  # size: (2M + 1) x K
    mtx_fri2amp_ri = linalg.lstsq(mtx_amp2freq_ri, np.eye(L))[0]
    # projection mtx_freq2visi to the null space of mtx_fri2amp
    mtx_null_proj = np.eye(L) - np.dot(mtx_fri2amp_ri.T,
                                       linalg.lstsq(mtx_fri2amp_ri.T, np.eye(L))[0])
    G_updated = np.dot(mtx_amp2visi_ri,
                       linalg.block_diag(*([mtx_fri2amp_ri] * num_bands))
                       ) + \
                np.dot(mtx_fri2visi_ri,
                       linalg.block_diag(*([mtx_null_proj] * num_bands))
                       )
    return G_updated 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def solve(self, Y):
    X, _residues, _rank, _sv = la.lstsq(self._Sigma, Y)
    return X 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def __truediv__(self, y, niter=100):
        if self.explicit is True:
            if isinstance(self.A, np.ndarray):
                if self.A.shape[0] == self.A.shape[1]:
                    xest = solve(self.A, y)
                else:
                    xest = lstsq(self.A, y)[0]
            else:
                xest = spsolve(self.A, y)
        else:
            xest = lsqr(self, y, iter_lim=niter)[0]
        return xest 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def apply(self, resp_matrix, orbit, weights=None):
        if weights is None:
            weights = np.eye(len(orbit))
        resp_matrix = np.dot(weights, resp_matrix)
        misallign = np.dot(weights, orbit)

        U, s, V = svd(resp_matrix)
        s_inv = np.zeros(len(s))
        s_max = max(s)
        for i in range(len(s)):
            #print("S[",i,"]=", s[i], "s max = ", s_max)
            if i < int(len(s)/2.):
                epsilon = self.epsilon_x
            else:
                epsilon = self.epsilon_y
            if s[i] <= s_max * epsilon:
                s_inv[i] = 0.
            else:
                s_inv[i] = 1. / s[i]
        Sinv = np.zeros((np.shape(U)[0], np.shape(V)[0]))
        Sinv[:len(s), :len(s)] = np.diag(s_inv)
        Sinv = np.transpose(Sinv)
        A = np.dot(np.transpose(V), np.dot(Sinv, np.transpose(U)))
        angle = np.dot(A, misallign)
        #a, b, _, _ = np.linalg.lstsq(resp_matrix, misallign, rcond=self.epsilon_x)
        #print(np.shape(a), np.shape(angle))
        #print(f"angle = {angle}, a = {a}")
        logger.debug("max(abs(angle)) = " + str(np.max(np.abs(angle))) + " min(abs(angle)) = " + str(np.min(np.abs(angle))))
        return angle 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def solver_lstsq(self, resp_matrix, orbit, weights):
        if weights is not None:
            # weights = np.eye(len(orbit))
            resp_matrix = np.dot(weights, resp_matrix)
            orbit = np.copy(np.dot(weights, orbit))
            
        #misallign = np.dot(weights, orbit)
        if np.shape(resp_matrix)[1] == 1:
            A = np.dot(resp_matrix.T, resp_matrix)
            Ainv = 1./A if A[0, 0] != 0 else A
            a = np.dot(np.dot(Ainv, resp_matrix.T), orbit)
        else:
            a, _, _, _ = np.linalg.lstsq(resp_matrix, orbit, rcond=self.epsilon_x)
            #a, _, _, _ = lstsq(resp_matrix, misallign, cond=self.epsilon_x, lapack_driver="gelsy")
        return a 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def glm(data, para):
    """
    Calculates a linear model fit ...
    anova/ancova/lin-regress/t-test/etc. Taken from:

    Returns
    -------
    statistic, p-value ???

    References
    ----------
    Peterson et al. Statistical limitations in functional neuroimaging
    I. Non-inferential methods and statistical models.  Phil Trans Royal Soc
    Lond B 354: 1239-1260.

    """
    if len(para) != len(data):
        raise ValueError("data and para must be same length in aglm")
    n = len(para)
    p = _support.unique(para)
    x = zeros((n,len(p)))  # design matrix
    for l in range(len(p)):
        x[:,l] = para == p[l]
    # fixme: normal equations are bad. Use linalg.lstsq instead.
    b = dot(dot(linalg.inv(dot(np.transpose(x),x)),  # i.e., b=inv(X'X)X'Y
                    np.transpose(x)),data)
    diffs = (data - dot(x,b))
    s_sq = 1./(n-len(p)) * dot(np.transpose(diffs), diffs)

    if len(p) == 2:  # ttest_ind
        c = array([1,-1])
        df = n-2
        fact = np.sum(1.0/np.sum(x,0),axis=0)  # i.e., 1/n1 + 1/n2 + 1/n3 ...
        t = dot(c,b) / np.sqrt(s_sq*fact)
        probs = betai(0.5*df,0.5,float(df)/(df+t*t))
        return t, probs
    else:
        raise ValueError("only ttest_ind implemented") 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def _solve_lsqr(self, X, y, shrinkage):
        """Least squares solver.

        The least squares solver computes a straightforward solution of the
        optimal decision rule based directly on the discriminant functions. It
        can only be used for classification (with optional shrinkage), because
        estimation of eigenvectors is not performed. Therefore, dimensionality
        reduction with the transform is not supported.

        Parameters
        ----------
        X : array-like, shape (n_samples, n_features)
            Training data.

        y : array-like, shape (n_samples,) or (n_samples, n_classes)
            Target values.

        shrinkage : string or float, optional
            Shrinkage parameter, possible values:
              - None: no shrinkage (default).
              - 'auto': automatic shrinkage using the Ledoit-Wolf lemma.
              - float between 0 and 1: fixed shrinkage parameter.

        Notes
        -----
        This solver is based on [1]_, section 2.6.2, pp. 39-41.

        References
        ----------
        .. [1] R. O. Duda, P. E. Hart, D. G. Stork. Pattern Classification
           (Second Edition). John Wiley & Sons, Inc., New York, 2001. ISBN
           0-471-05669-3.
        """
        self.means_ = _class_means(X, y)
        self.covariance_ = _class_cov(X, y, self.priors_, shrinkage)
        self.coef_ = linalg.lstsq(self.covariance_, self.means_.T)[0].T
        self.intercept_ = (-0.5 * np.diag(np.dot(self.means_, self.coef_.T)) +
                           np.log(self.priors_)) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def test_multinomial_grad_hess():
    rng = np.random.RandomState(0)
    n_samples, n_features, n_classes = 100, 5, 3
    X = rng.randn(n_samples, n_features)
    w = rng.rand(n_classes, n_features)
    Y = np.zeros((n_samples, n_classes))
    ind = np.argmax(np.dot(X, w.T), axis=1)
    Y[range(0, n_samples), ind] = 1
    w = w.ravel()
    sample_weights = np.ones(X.shape[0])
    grad, hessp = _multinomial_grad_hess(w, X, Y, alpha=1.,
                                         sample_weight=sample_weights)
    # extract first column of hessian matrix
    vec = np.zeros(n_features * n_classes)
    vec[0] = 1
    hess_col = hessp(vec)

    # Estimate hessian using least squares as done in
    # test_logistic_grad_hess
    e = 1e-3
    d_x = np.linspace(-e, e, 30)
    d_grad = np.array([
        _multinomial_grad_hess(w + t * vec, X, Y, alpha=1.,
                               sample_weight=sample_weights)[0]
        for t in d_x
    ])
    d_grad -= d_grad.mean(axis=0)
    approx_hess_col = linalg.lstsq(d_x[:, np.newaxis], d_grad)[0].ravel()
    assert_array_almost_equal(hess_col, approx_hess_col) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def fit(self):
        """Solve linear system for the weights.

        The weights  `self.w` (:math:`\mathbf w`) are found from: :math:`\mathbf
        G\,\mathbf w = \mathbf z` or if :math:`r` is given :math:`(\mathbf G +
        r\,\mathbf I)\,\mathbf w = \mathbf z`.

        with centers == points (center vectors are all data points). Then G is
        quadratic. Updates ``self.w``.

        Notes
        -----
        ``self.r != None`` : linear system solver
            Use :func:`scipy.linalg.solve`. For :math:`r=0`, this always yields
            perfect interpolation at the data points. May be numerically
            unstable in that case. Use :math:`r>0` to increase stability (try
            small values such as ``1e-10`` first) or create smooth fitting (generate
            more stiff functions with higher `r`). Behaves similar to ``lstsq``
            but appears to be numerically more stable (no small noise in
            solution) .. but `r` it is another parameter that needs to be
            tuned.
        ``self.r = None`` : least squares solver
            Use :func:`scipy.linalg.lstsq`. Numerically more stable than
            direct solver w/o regularization. Will mostly be the same as
            the interpolation result, but will not go thru all points for
            very noisy data. May create small noise in solution (plot fit
            with high point density).
        """
        # this test may be expensive for big data sets
        assert (self.centers == self.points).all(), "centers == points not fulfilled"
        # re-use self.distsq if possible
        if self.distsq is None:
            self.distsq = self.get_distsq()
        G = self.rbf(self.distsq, self.p)
        if self.r is None:
            self.w, res, rnk, svs = linalg.lstsq(G, self.values)
        else:
            self.w = linalg.solve(G + np.identity(G.shape[0])*self.r, self.values) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def post(self, hyp, X, y):
        """ Generic function to compute posterior distribution.

            This function will save the posterior mean and precision matrix as
            self.m and self.A and will also update internal parameters (e.g.
            N, D and the prior covariance (Sigma_a) and precision (Lambda_a).
        """

        N = X.shape[0]
        if len(X.shape) == 1:
            D = 1
        else:
            D = X.shape[1]

        if (hyp == self.hyp).all() and hasattr(self, 'N'):
            print("hyperparameters have not changed, exiting")
            return
        
        beta, alpha, gamma = self._parse_hyps(hyp, X)

        if self.verbose:
            print("estimating posterior ... | hyp=", hyp)

        # prior variance
        if len(alpha) == 1 or len(alpha) == D:
            self.Sigma_a = np.diag(np.ones(D))/alpha
            self.Lambda_a = np.diag(np.ones(D))*alpha
        else:
            raise ValueError("hyperparameter vector has invalid length")

        # compute posterior precision and mean
        self.A = X.T.dot(self.Lambda_n).dot(X) + self.Lambda_a
        self.m = linalg.solve(self.A, X.T, 
                              check_finite=False).dot(self.Lambda_n).dot(y)
        #self.m = linalg.lstsq(self.A, X.T, 
        #                      check_finite=False)[0].dot(self.Lambda_n).dot(y)

        # save stuff
        self.N = N
        self.D = D
        self.hyp = hyp 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import lstsq [as 别名]
def fit_line(x, y, plot_line=False):
    """
    Plot best fit least squares line, return R^2
    """
    A = vstack([x, ones(len(x))]).T
    m, c = linalg.lstsq(A, y)[0]
    if plot_line:
        plt.plot(x, m*x + c, 'r', lw=1.2)
    return square(stats.pearsonr(x, y)), m, c