Python linalg.cho_factor方法代码示例

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def mvn_loglike(x, sigma):
    '''loglike multivariate normal

    assumes x is 1d, (nobs,) and sigma is 2d (nobs, nobs)

    brute force from formula
    no checking of correct inputs
    use of inv and log-det should be replace with something more efficient
    '''
    #see numpy thread
    #Sturla: sqmahal = (cx*cho_solve(cho_factor(S),cx.T).T).sum(axis=1)
    sigmainv = linalg.inv(sigma)
    logdetsigma = np.log(np.linalg.det(sigma))
    nobs = len(x)

    llf = - np.dot(x, np.dot(sigmainv, x))
    llf -= nobs * np.log(2 * np.pi)
    llf -= logdetsigma
    llf *= 0.5
    return llf 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def do_sampleF2(Y, X, current_F1, log_theta):

    log_theta_sigma2, log_theta_lengthscale, log_theta_lambda = unpack_log_theta(log_theta)

    n = Y.shape[0]

    K_F2 = covariance_function(current_F1, current_F1, log_theta_sigma2[1], log_theta_lengthscale[1]) + np.eye(n) * 1e-9
    K_Y = K_F2 + np.eye(n) * np.exp(log_theta_lambda)
    L_K_Y, lower_K_Y = linalg.cho_factor(K_Y, lower=True)
    K_inv_Y = linalg.cho_solve((L_K_Y, lower_K_Y), Y)

    mu = np.dot(K_F2, K_inv_Y)

    K_inv_K = linalg.cho_solve((L_K_Y, lower_K_Y), K_F2)
    Sigma = K_F2 - np.dot(K_F2, K_inv_K)

    L_Sigma = linalg.cholesky(Sigma, lower=True)
    proposed_F2 = mu + np.dot(L_Sigma, np.random.normal(0.0, 1.0, [n, 1]))

    return proposed_F2

## Unpack the vector of parameters into their three elements 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def inv(self, logdet=False):
        if self.ndim == 1:
            inv = 1.0 / self

            if logdet:
                return inv, np.sum(np.log(self))
            else:
                return inv
        else:
            try:
                cf = sl.cho_factor(self)
                inv = sl.cho_solve(cf, np.identity(cf[0].shape[0]))
                if logdet:
                    ld = 2.0 * np.sum(np.log(np.diag(cf[0])))
            except np.linalg.LinAlgError:
                u, s, v = np.linalg.svd(self)
                inv = np.dot(u / s, u.T)
                if logdet:
                    ld = np.sum(np.log(s))
            if logdet:
                return inv, ld
            else:
                return inv 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def standardized_forecasts_error(self):
        """
        Standardized forecast errors
        """
        if self._standardized_forecasts_error is None:
            from scipy import linalg
            self._standardized_forecasts_error = np.zeros(
                self.forecasts_error.shape, dtype=self.dtype)

            for t in range(self.forecasts_error_cov.shape[2]):
                if self.nmissing[t] > 0:
                    self._standardized_forecasts_error[:, t] = np.nan
                if self.nmissing[t] < self.k_endog:
                    mask = ~self.missing[:, t].astype(bool)
                    F = self.forecasts_error_cov[np.ix_(mask, mask, [t])]
                    upper, _ = linalg.cho_factor(F[:, :, 0])
                    self._standardized_forecasts_error[mask, t] = (
                        linalg.solve_triangular(
                            upper, self.forecasts_error[mask, t]
                        )
                    )

        return self._standardized_forecasts_error 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def log_likelihood(self, p):
        p = self.to_params(p)

        v = self.rvs(p)

        res = self.vs - v - p['mu']

        cov = p['nu']*p['nu']*np.diag(self.dvs*self.dvs)
        cov += generate_covariance(self.ts, p['sigma'], p['tau'])

        cfactor = sl.cho_factor(cov)
        cc, lower = cfactor

        n = self.ts.shape[0]

        return -0.5*n*np.log(2.0*np.pi) - np.sum(np.log(np.diag(cc))) - 0.5*np.dot(res, sl.cho_solve(cfactor, res)) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def __init__(self, data):
        self._data = np.atleast_2d(data)

        self._mean = np.mean(data, axis=0)
        self._cov = None

        if self.data.shape[0] > 1:
            try:
                self._cov = np.cov(data.T)

                # Try factoring now to see if regularization is needed
                la.cho_factor(self._cov)

            except la.LinAlgError:
                self._cov = oas_cov(data)

        self._set_bandwidth()

        # store transformation variables for drawing random values
        alphas = np.std(data, axis=0)
        ms = 1./alphas
        m_i, m_j = np.meshgrid(ms, ms)
        ms = m_i * m_j
        self._draw_cov = ms * self._kernel_cov
        self._scale_fac = alphas 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def _set_bandwidth(self):
        r"""
        Use Scott's rule to set the kernel bandwidth:

        .. math::

            \mathcal{K} = n^{-1/(d+4)} \Sigma^{1/2}

        Also store Cholesky decomposition for later.
        """
        if self.size > 0 and self._cov is not None:
            self._kernel_cov = self._cov * self.size ** (-2/(self.ndim + 4))

            # Used to evaluate PDF with cho_solve()
            self._cho_factor = la.cho_factor(self._kernel_cov)

            # Make sure the estimated PDF integrates to 1.0
            self._lognorm = self.ndim/2 * np.log(2*np.pi) + np.log(self.size) +\
                np.sum(np.log(np.diag(self._cho_factor[0])))

        else:
            self._lognorm = -np.inf 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def _evaluate_point_logpdf(args):
    """
    Evaluate the Gaussian KDE at a given point `p`.  This lives outside the KDE method to allow for
    parallelization using :mod:`multipocessing`. Since :func:`map` only allows single-argument
    functions, the following arguments to be packed into a single tuple.

    :param p:
        The point to evaluate the KDE at.

    :param data:
        The `(N, ndim)`-shaped array of data used to construct the KDE.

    :param cho_factor:
        A Cholesky decomposition of the kernel covariance matrix.
    """
    point, data, cho_factor = args

    # Use Cholesky decomposition to avoid direct inversion of covariance matrix
    diff = data - point
    tdiff = la.cho_solve(cho_factor, diff.T, check_finite=False).T
    diff *= tdiff

    # Work in the log to avoid large numbers
    return logsumexp(-np.sum(diff, axis=1)/2) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def __init__(self, Sigma):
    self._cholesky = la.cho_factor(Sigma) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def mvn_nloglike_obs(x, sigma):
    '''loglike multivariate normal

    assumes x is 1d, (nobs,) and sigma is 2d (nobs, nobs)

    brute force from formula
    no checking of correct inputs
    use of inv and log-det should be replace with something more efficient
    '''
    #see numpy thread
    #Sturla: sqmahal = (cx*cho_solve(cho_factor(S),cx.T).T).sum(axis=1)

    #Still wasteful to calculate pinv first
    sigmainv = linalg.inv(sigma)
    cholsigmainv = linalg.cholesky(sigmainv)
    #2 * np.sum(np.log(np.diagonal(np.linalg.cholesky(A)))) #Dag mailinglist
    # logdet not needed ???
    #logdetsigma = 2 * np.sum(np.log(np.diagonal(cholsigmainv)))
    x_whitened = np.dot(cholsigmainv, x)

    #sigmainv = linalg.cholesky(sigma)
    logdetsigma = np.log(np.linalg.det(sigma))

    sigma2 = 1. # error variance is included in sigma

    llike  =  0.5 * (np.log(sigma2) - 2.* np.log(np.diagonal(cholsigmainv))
                          + (x_whitened**2)/sigma2
                          +  np.log(2*np.pi))

    return llike, (x_whitened**2) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def _make_AB_matrices(self):
		# Vertical
		n = self.num_stencils_vertical
		cov_zz = self.cov_matrix_vertical[:n,:n]
		cov_xz = self.cov_matrix_vertical[n:, :n]
		cov_zx = self.cov_matrix_vertical[:n, n:]
		cov_xx = self.cov_matrix_vertical[n:, n:]

		cf = linalg.cho_factor(cov_zz)
		inv_cov_zz = linalg.cho_solve(cf, np.eye(cov_zz.shape[0]))

		self.A_vertical = cov_xz.dot(inv_cov_zz)

		BBt = cov_xx - self.A_vertical.dot(cov_zx)

		U, S, Vt = np.linalg.svd(BBt)
		L = np.sqrt(S[:self.input_grid.dims[0]])

		self.B_vertical = U * L

		# Horizontal
		n = self.num_stencils_horizontal
		cov_zz = self.cov_matrix_horizontal[:n,:n]
		cov_xz = self.cov_matrix_horizontal[n:, :n]
		cov_zx = self.cov_matrix_horizontal[:n, n:]
		cov_xx = self.cov_matrix_horizontal[n:, n:]

		cf = linalg.cho_factor(cov_zz)
		inv_cov_zz = linalg.cho_solve(cf, np.eye(cov_zz.shape[0]))

		self.A_horizontal = cov_xz.dot(inv_cov_zz)

		BBt = cov_xx - self.A_horizontal.dot(cov_zx)

		U, S, Vt = np.linalg.svd(BBt)
		L = np.sqrt(S[:self.input_grid.dims[1]])

		self.B_horizontal = U * L 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def makeAMatrix(self):
        """
        Calculates the "A" matrix, that uses the existing data to find a new 
        component of the new phase vector.
        """
        # Cholsky solve can fail - if so do brute force inversion
        try:
            cf = linalg.cho_factor(self.cov_mat_zz)
            inv_cov_zz = linalg.cho_solve(cf, numpy.identity(self.cov_mat_zz.shape[0]))
        except linalg.LinAlgError:
            raise linalg.LinAlgError("Could not invert Covariance Matrix to for A and B Matrices. Try with a larger pixel scale")

        self.A_mat = self.cov_mat_xz.dot(inv_cov_zz) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def log_p_Y_given_F1(Y, F1, log_theta):

    log_theta_sigma2, log_theta_lengthscale, log_theta_lambda = unpack_log_theta(log_theta)

    n = Y.shape[0]

    K_Y = covariance_function(F1, F1, log_theta_sigma2[1], log_theta_lengthscale[1]) + np.eye(n) * np.exp(log_theta_lambda)
    L_K_Y, lower_K_Y = linalg.cho_factor(K_Y, lower=True)

    nu = linalg.solve_triangular(L_K_Y, Y, lower=True)

    return -np.sum(np.log(np.diagonal(L_K_Y))) - 0.5 * np.dot(nu.transpose(), nu)


## Elliptical Slice Sampling to sample from the posterior over latent variables at layer 1 (the latent variables at layer 2 are integrated out analytically) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def rake_mvdr_filters(self, source, interferer, R_n, delay=0.03, epsilon=5e-3):
        """
        Compute the time-domain filters of the minimum variance distortionless
        response beamformer.
        """

        H = build_rir_matrix(
            self.R,
            (source, interferer),
            self.Lg,
            self.fs,
            epsilon=epsilon,
            unit_damping=True,
        )
        L = H.shape[1] // 2

        # the constraint vector
        kappa = int(delay * self.fs)
        h = H[:, kappa]

        # We first assume the sample are uncorrelated
        R_xx = np.dot(H[:, :L], H[:, :L].T)
        K_nq = np.dot(H[:, L:], H[:, L:].T) + R_n

        # Compute the TD filters
        C = la.cho_factor(R_xx + K_nq, check_finite=False)
        g_val = la.cho_solve(C, h)

        g_val /= np.inner(h, g_val)
        self.filters = g_val.reshape((self.M, self.Lg))

        # compute and return SNR
        num = np.inner(g_val.T, np.dot(R_xx, g_val))
        denom = np.inner(np.dot(g_val.T, K_nq), g_val)

        return num / denom 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import cho_factor [as 别名]
def solve(self, other):
        if other.ndim == 1:
            Nx = np.array(other / self.N)
        elif other.ndim == 2:
            Nx = np.array(other / self.N[:, None])
        UNx = np.dot(self.U.T, Nx)

        Sigma = np.diag(1 / self.J) + np.dot(self.U.T, self.U / self.N[:, None])
        cf = sl.cho_factor(Sigma)
        if UNx.ndim == 1:
            tmp = np.dot(self.U, sl.cho_solve(cf, UNx)) / self.N
        else:
            tmp = np.dot(self.U, sl.cho_solve(cf, UNx)) / self.N[:, None]
        return Nx - tmp