我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.linalg.inv()。
def __call__(self, params): print '???', params sd1 = params[0] sd2 = params[1] cor = params[2] if sd1 < 0. or sd1 > 10. or sd2 < 0. or sd2 > 10. or cor < -1. or cor > 1.: return np.inf bandwidth = maths.stats.choleskysqrt2d(sd1, sd2, cor) bandwidthdet = la.det(bandwidth) bandwidthinv = la.inv(bandwidth) diff = sample[self.__iidx] - sample[self.__jidx] temp = diff.dot(bandwidthinv.T) temp *= temp e = np.exp(np.sum(temp, axis=1)) s = np.sum(e**(-.25) - 4 * e**(-.5)) cost = self.__n / bandwidthdet + (2. / bandwidthdet) * s print '!!!', cost return cost / 10000.
def distribution_parameters(self, parameter_name): if parameter_name=='trend': E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix) mean = dot(inv(eye(self.size)+E),self.data) cov = (self.parameters.list['sigma2'].current_value)*inv(eye(self.size)+E) return {'mean' : mean, 'cov' : cov} elif parameter_name=='sigma2': E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix) pos = self.size loc = 0 scale = 0.5*dot((self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)).T,(self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)))+0.5*dot(dot(self.parameters.list['trend'].current_value.T,E),self.parameters.list['trend'].current_value) elif parameter_name=='lambda2': pos = self.size-self.total_variation_order-1+self.alpha loc = 0.5*(norm(dot(self.derivative_matrix,self.parameters.list['trend'].current_value),ord=1))/self.parameters.list['sigma2'].current_value+self.rho scale = 1 elif parameter_name==str('omega'): pos = [sqrt(((self.parameters.list['lambda2'].current_value**2)*self.parameters.list['sigma2'].current_value)/(dj**2)) for dj in dot(self.derivative_matrix,self.parameters.list['trend'].current_value)] loc = 0 scale = self.parameters.list['lambda2'].current_value**2 return {'pos' : pos, 'loc' : loc, 'scale' : scale}
def test_byteorder_check(): # Byte order check should pass for native order if sys.byteorder == 'little': native = '<' else: native = '>' for dtt in (np.float32, np.float64): arr = np.eye(4, dtype=dtt) n_arr = arr.newbyteorder(native) sw_arr = arr.newbyteorder('S').byteswap() assert_equal(arr.dtype.byteorder, '=') for routine in (linalg.inv, linalg.det, linalg.pinv): # Normal call res = routine(arr) # Native but not '=' assert_array_equal(res, routine(n_arr)) # Swapped assert_array_equal(res, routine(sw_arr))
def test_basic(self): import numpy.linalg as linalg A = np.array([[1., 2.], [3., 4.]]) mA = matrix(A) assert_(np.allclose(linalg.inv(A), mA.I)) assert_(np.all(np.array(np.transpose(A) == mA.T))) assert_(np.all(np.array(np.transpose(A) == mA.H))) assert_(np.all(A == mA.A)) B = A + 2j*A mB = matrix(B) assert_(np.allclose(linalg.inv(B), mB.I)) assert_(np.all(np.array(np.transpose(B) == mB.T))) assert_(np.all(np.array(np.transpose(B).conj() == mB.H)))
def test_basic(self): import numpy.linalg as linalg A = np.array([[1., 2.], [3., 4.]]) mA = matrix(A) B = np.identity(2) for i in range(6): assert_(np.allclose((mA ** i).A, B)) B = np.dot(B, A) Ainv = linalg.inv(A) B = np.identity(2) for i in range(6): assert_(np.allclose((mA ** -i).A, B)) B = np.dot(B, Ainv) assert_(np.allclose((mA * mA).A, np.dot(A, A))) assert_(np.allclose((mA + mA).A, (A + A))) assert_(np.allclose((3*mA).A, (3*A))) mA2 = matrix(A) mA2 *= 3 assert_(np.allclose(mA2.A, 3*A))
def test_model_pq(self): """Test the model PQ matrix generation """ new_model = self.bayes.update(models={ 'simp': self.model.update_dof([4, 2])}) P, q = new_model._get_model_pq() epsilon = np.array([2, 1]) sigma = inv(np.diag(np.ones(2))) P_true = sigma q_true = -np.dot(epsilon, sigma) npt.assert_array_almost_equal(P, P_true, decimal=8) npt.assert_array_almost_equal(q, q_true, decimal=8)
def test_sim_pq(self): """Test the simulation PQ matrix generation """ initial_data = self.bayes.get_data() sens_matrix = self.bayes._get_sens(initial_data=initial_data) P, q = self.bayes._get_sim_pq(initial_data, sens_matrix) sigma = self.exp1.get_sigma() P_true = np.dot(np.dot(sens_matrix['simple'].T, inv(sigma)), sens_matrix['simple']) npt.assert_array_almost_equal(P, P_true, decimal=8) epsilon = self.bayes.simulations['simple']['exp']\ .compare(initial_data['simple']) q_true = -np.dot(np.dot(epsilon, inv(sigma)), sens_matrix['simple']) npt.assert_array_almost_equal(q, q_true, decimal=8)
def get_log_like(self, sim_data): """Gets the log likelihood of the current simulation Args: sim_data(list): Length three list corresponding to the `__call__` from a Experiment object experiment(Experiment): A valid Experiment object Return: (float): The log of the likelihood of the simulation """ epsilon = self.compare(sim_data) return -0.5 * np.dot(epsilon, np.dot(inv(self.get_sigma()), epsilon))
def get_pq(self, scale=False): """Returns the P and q matrix components for the model Keyword Args: scale(bool): Flag to scale the model values """ prior_scale = self.get_scaling() prior_var = inv(self.get_sigma()) prior_delta = self.get_dof() - self.prior.get_dof() if scale: return np.dot(prior_scale, np.dot(prior_var, prior_scale)),\ -np.dot(prior_scale, np.dot(prior_delta, prior_var)) else: return prior_var, -np.dot(prior_delta, prior_var) # end
def trainClassifer(self,labels,vectors,verbose=False,ilog=None): ''' Do not call this function instead call train. ''' self.training_size = len(labels) c = len(labels) r = len(vectors[0]) y = array(labels,'d') X = zeros((r,c),'d') for i in range(len(vectors)): X[:,i] = vectors[i] tmp1 = inv(self.lam*eye(r) + dot(X,X.transpose())) tmp2 = dot(y,X.transpose()) self.w = w = dot(tmp1,tmp2)
def trainRidgeRegression(self,labels,vectors,verbose=False): self.training_size = len(labels) c = len(labels) r = len(vectors[0]) self.c = c self.r = r y = array(labels,'d') X = zeros((r,c),'d') for i in range(len(vectors)): X[:,i] = vectors[i] self.X = X kernel_matrix = zeros((c,c),'d') for i in range(c): kernel_matrix[:,i] = self.kernel(X,X[:,i:i+1]) self.w = w = dot(y,inv(kernel_matrix + self.lam*eye(c)))
def coefficient_variance(number_of_coefficient, y_arr, y_explained, x_matrix): """ :param number_of_coefficient: number from range(n) :param y_arr: input y :param y_explained: y explained array :param x_matrix: x matrix from docs :return: variance of coefficient """ variance_e_2 = rss(y_arr, y_explained) / (n - k) # (X^T * X)^-1 mmatrix = linal.inv(numpy.dot(x_matrix.T, x_matrix)) v_matrix = numpy.dot(variance_e_2, mmatrix) return v_matrix[number_of_coefficient][number_of_coefficient] ############################################################## ######### ??????? ???????? ??????????? y ?? x2, x3 ? x4####### ####### ??????? ?????????? ????????? ######################### ##############################################################
def fast_optimize(endog, exog, n_obs=0, n_vars=0, max_iter=10000, tolerance=1e-10): """ A convenience function for the Newton-Raphson method to evaluate a logistic model. :param endog: An Nx1 vector of endogenous predictions :param exog: An NxK vector of exogenous predictors :param n_obs: The number of observations N :param n_vars: The number of exogenous predictors K :param max_iter: The maximum number of iterations :param tolerance: Margin of error for convergence :return: The error-minimizing parameters for the model. """ iterations = 0 oldparams = np.inf newparams = np.repeat(0, n_vars) while iterations < max_iter and np.any(np.abs(newparams - oldparams) > tolerance): oldparams = newparams try: H = logit_hessian(exog, oldparams, n_obs) newparams = oldparams - dot(inv(H), logit_score(endog, exog, oldparams, n_obs)) except LinAlgError: raise LinAlgError iterations += 1 return newparams
def update(self, z): sigmas_f, sigmas_h = self.sigmas_f, self.sigmas_h for i in range(self.M.num_sigmas): sigmas_h[i] = self.h(sigmas_f[i]) zp, Pz = self.unscented_transform(sigmas_h, self.M.Wm, self.M.Wc, self._R) Pxz = np.zeros((self._n, self._k)) for i in range(self.M.num_sigmas): Pxz += self.M.Wc[i] * np.outer(sigmas_f[i] - self.xp, sigmas_h[i] - zp) K = Pxz.dot(inv(Pz)) # Kalman gain self._x = self.xp + K.dot(z-zp) self._P = self.Pp - K.dot(Pz).dot(K.T)
def get_dimensions(self): # Calculate within class scatter Sw = np.sum(self._scatter_matrix_by_class.values(), axis=0) # Calculate between class scatter s = (self._p, self._p) Sb = np.zeros(s) for k in self._classes: a = self._mean_vector_by_class[k] - self._global_mean_vector Sb += self._N_by_class[k] * a.dot(a.T) # Compute eigenvectors Sw_inv = inv(Sw) A = Sw_inv.dot(Sb) eigen_values, eigen_vectors = eig(A) idx = np.argsort(eigen_values) eigen_vectors = eigen_vectors[idx][::-1] return eigen_vectors
def _udpate_item_features(self): # Gibbs sampling for item features for item_id in xrange(self.n_item): indices = self.ratings_csc_[:, item_id].indices features = self.user_features_[indices, :] rating = self.ratings_csc_[:, item_id].data - self.mean_rating_ rating = np.reshape(rating, (rating.shape[0], 1)) covar = inv(self.alpha_item + self.beta * np.dot(features.T, features)) lam = cholesky(covar) temp = (self.beta * np.dot(features.T, rating) + np.dot(self.alpha_item, self.mu_item)) mean = np.dot(covar, temp) temp_feature = mean + np.dot( lam, self.rand_state.randn(self.n_feature, 1)) self.item_features_[item_id, :] = temp_feature.ravel()
def _update_user_features(self): # Gibbs sampling for user features for user_id in xrange(self.n_user): indices = self.ratings_csr_[user_id, :].indices features = self.item_features_[indices, :] rating = self.ratings_csr_[user_id, :].data - self.mean_rating_ rating = np.reshape(rating, (rating.shape[0], 1)) covar = inv( self.alpha_user + self.beta * np.dot(features.T, features)) lam = cholesky(covar) # aplha * sum(V_j * R_ij) + LAMBDA_U * mu_u temp = (self.beta * np.dot(features.T, rating) + np.dot(self.alpha_user, self.mu_user)) # mu_i_star mean = np.dot(covar, temp) temp_feature = mean + np.dot( lam, self.rand_state.randn(self.n_feature, 1)) self.user_features_[user_id, :] = temp_feature.ravel()
def _update_user_feature(self): """Fix item features and update user features """ for i in xrange(self.n_user): _, item_idx = self.ratings_csr_[i, :].nonzero() # number of ratings of user i n_u = item_idx.shape[0] if n_u == 0: logger.debug("no ratings for user %d", i) continue item_features = self.item_features_.take(item_idx, axis=0) ratings = self.ratings_csr_[i, :].data - self.mean_rating_ A_i = (np.dot(item_features.T, item_features) + self.reg * n_u * np.eye(self.n_feature)) V_i = np.dot(item_features.T, ratings) self.user_features_[i, :] = np.dot(inv(A_i), V_i)
def warpImage(background, image, parallelogram): ''' @params: background is unchanged image image is image to be warped parallelogram is the coordinates to warp the image to, starting at upper left and going clockwise returns a new image that is the composition of background and image after image has been warped ''' mapped = np.array([[parallelogram[0].x, parallelogram[1].x, parallelogram[2].x], [parallelogram[0].y, parallelogram[1].y, parallelogram[2].y], [1,1,1]]) width, height = image.size original = np.array([[0, width, width],[0, 0, height]]) #solve for affine matrix solution = np.dot(original, inv(mapped)) #unroll matrix into a sequence affine = (solution[0][0], solution[0][1], solution[0][2], solution[1][0], solution[1][1], solution[1][2]) transformed = image.transform(background.size, Image.AFFINE, affine) white = Image.new("RGBA", (width, height), "white") transformedMask = white.transform(background.size, Image.AFFINE, affine) background.paste(transformed, (0,0), transformedMask) return background
def tforminv(trans, uv): """ Function: ---------- apply the inverse of affine transform 'trans' to uv Parameters: ---------- @trans: 3x3 np.array transform matrix @uv: Kx2 np.array each row is a pair of coordinates (x, y) Returns: ---------- @xy: Kx2 np.array each row is a pair of inverse-transformed coordinates (x, y) """ Tinv = inv(trans) xy = tformfwd(Tinv, uv) return xy
def aryule(c, k): """Solve Yule-Walker equation. Args: c (numpy array): Coefficients (i.e. autocorrelation) k (int): Assuming the AR(k) model Returns: numpy array: k model parameters Some formulations solve: C a = -c, but we actually solve C a = c. """ a = np.zeros(k) # ignore a singular matrix C = toeplitz(c[:k]) if not np.all(C == 0.0) and np.isfinite(ln.cond(C)): a = np.dot(ln.inv(C), c[1:]) return a
def estimate_X_TLS(y,H): """ Estimator of x for the Linear Model using Total Least Square (TLS). As compared to the classical Least Squares Estimator, the TLS estimator is more suited when H is not exactly known or contained with noise [GOL80]_. .. [GOL80] Golub, Gene H., and Charles F. Van Loan. "An analysis of the total least squares problem." SIAM Journal on Numerical Analysis 17.6 (1980): 883-893. """ N,L=H.shape H=np.matrix(H) Z=np.hstack([H,y]) U,S,V=lg.svd(Z) V=np.matrix(V) V=V.H VHY=V[:L,L:]; VYY=V[L:,L:]; x= -VHY*lg.inv(VYY); return x
def make_eigvals_positive(am, targetprod): """For the symmetric square matrix `am`, increase any zero eigenvalues such that the total product of eigenvalues is greater or equal to `targetprod`. Returns a (possibly) new, non-singular matrix.""" w, v = linalg.eigh(am) # use eigh since a is symmetric mask = w < 1.e-10 if np.any(mask): nzprod = np.product(w[~mask]) # product of nonzero eigenvalues nzeros = mask.sum() # number of zero eigenvalues new_val = max(1.e-10, (targetprod / nzprod) ** (1. / nzeros)) w[mask] = new_val # adjust zero eigvals am_new = np.dot(np.dot(v, np.diag(w)), linalg.inv(v)) # re-form cov else: am_new = am return am_new
def compSpecSample(self, angle): """ Computes the spectrum in a sample with Capon beamformer. MUST compute the covariance matrix (call `compCov()`) before calling this function. Parameters ---------- angle : float Direction-of-arrival (DOA) angle in range [0, pi). Returns ------- p : float The response of Capon beamformer. """ a = self.sarr.steer(angle) p = 1. / (a.T.conj().dot(inv(self.r).dot(a))) # Discard imaginary part return p.real
def compSpecSample(self, angle): """ Computes the spatial spectrum response in a specified angle with FLOM matrix. MUST compute the FLOM matrix (call `compFLOM()`) before calling this function. Parameters ---------- angle : float Direction-of-arrival (DOA) angle in range [0, pi). Returns ------- p : float The response of spatial spectrum. """ a = self.sarr.steer(angle) p = 1. / (a.T.conj().dot(inv(self.gamma).dot(a))) # Discard imaginary part return p.real
def _mvee(self, points, tolerance): # Taken from: http://stackoverflow.com/questions/14016898/port-matlab-bounding-ellipsoid-code-to-python points = np.asmatrix(points) n, d = points.shape q = np.column_stack((points, np.ones(n))).T err = tolerance + 1.0 u = np.ones(n) / n while err > tolerance: # assert u.sum() == 1 # invariant x = q * np.diag(u) * q.T m = np.diag(q.T * la.inv(x) * q) jdx = np.argmax(m) step_size = (m[jdx] - d - 1.0) / ((d + 1) * (m[jdx] - 1.0)) new_u = (1 - step_size) * u new_u[jdx] += step_size err = la.norm(new_u - u) u = new_u c = u * points a = la.inv(points.T * np.diag(u) * points - c.T * c) / d return np.asarray(a), np.squeeze(np.asarray(c))
def _int_var_rbf_hyp(self, hyp, X, jitter=1e-8): """ Posterior integral variance as a function of hyper-parameters :param hyp: RBF kernel hyper-parameters [s2, el_1, ..., el_d] :param X: sigma-points :param jitter: numerical jitter (for stabilizing computations) :return: posterior integral variance """ # reshape X to SP matrix X = np.reshape(X, (self.n, self.d)) # set kernel hyper-parameters s2, el = 1, hyp # sig_var hyper always set to 1 self.kern.param_array[0] = s2 # variance self.kern.param_array[1:] = el # lengthscale K = self.kern.K(X) L = np.diag(el ** 2) # posterior variance of the integral ks = s2 * np.sqrt(det(L + np.eye(self.d))) * multivariate_normal(mean=np.zeros(self.d), cov=L).pdf(X) postvar = s2 * np.sqrt(det(2 * inv(L) + np.eye(self.d))) ** -1 - ks.dot( solve(K + jitter * np.eye(self.n), ks.T)) return postvar
def cov(self): """ Compute parameter covariance Returns ------- c : ndarray Parameter covariance """ s = self.s nobs = self._xe.shape[0] scale = 1 / (nobs - int(self._debiased) * self._df) if self.square: ji = self.inv_jacobian out = ji @ s @ ji.T else: j = self.jacobian out = inv(j.T @ inv(s) @ j) out = (scale / 2) * (out + out.T) return out
def w(self, moments): """ Score/moment condition weighting matrix Parameters ---------- moments : ndarray Moment conditions (nobs by nmoments) Returns ------- w : ndarray Weighting matrix computed from moment conditions """ if self._center: moments = moments - moments.mean(0)[None, :] nobs = moments.shape[0] out = moments.T @ moments / nobs return inv((out + out.T) / 2.0)
def w(self, moments): """ Score/moment condition weighting matrix Parameters ---------- moments : ndarray Moment conditions (nobs by nmoments) Returns ------- w : ndarray Weighting matrix computed from moment conditions """ if self._center: moments = moments - moments.mean(0)[None, :] out = self._kernel_cov(moments) return inv(out)
def _f_statistic(self, params, cov, debiased): non_const = ~(self._x.ptp(0) == 0) test_params = params[non_const] test_cov = cov[non_const][:, non_const] test_stat = test_params.T @ inv(test_cov) @ test_params test_stat = float(test_stat) nobs, nvar = self._x.shape null = 'All parameters ex. constant are zero' name = 'Model F-statistic' df = test_params.shape[0] if debiased: wald = WaldTestStatistic(test_stat / df, null, df, nobs - nvar, name=name) else: wald = WaldTestStatistic(test_stat, null, df, name=name) return wald
def cov(self): """Covariance of estimated parameters""" x, z = self.x, self.z nobs, nvar = x.shape pinvz = self._pinvz v = (x.T @ z) @ (pinvz @ x) / nobs if self._kappa != 1: kappa = self._kappa xpx = x.T @ x / nobs v = (1 - kappa) * xpx + kappa * v vinv = inv(v) c = vinv @ self.s @ vinv / nobs return (c + c.T) / 2
def _gls_cov(self): x = self._x sigma = self._sigma sigma_inv = inv(sigma) xpx = blocked_inner_prod(x, sigma_inv) # Handles case where sigma_inv is not inverse of full_sigma xeex = blocked_inner_prod(x, sigma_inv @ self._full_sigma @ sigma_inv) if self._constraints is None: xpxi = inv(xpx) cov = xpxi @ xeex @ xpxi else: cons = self._constraints xpx = cons.t.T @ xpx @ cons.t xpxi = inv(xpx) xeex = cons.t.T @ xeex @ cons.t cov = cons.t @ (xpxi @ xeex @ xpxi) @ cons.t.T cov = (cov + cov.T) / 2 return cov
def __init__(self, x, eps, sigma, full_sigma, gls=False, debiased=False, constraints=None): super(HeteroskedasticCovariance, self).__init__(x, eps, sigma, full_sigma, gls=gls, debiased=debiased, constraints=constraints) self._name = 'Heteroskedastic (Robust) Covariance' k = len(x) weights = inv(sigma) if gls else eye(k) bigx = blocked_diag_product(x, weights) nobs = eps.shape[0] e = eps.T.ravel()[:, None] bigxe = bigx * e m = bigx.shape[1] xe = zeros((nobs, m)) for i in range(nobs): xe[i, :] = bigxe[i::nobs].sum(0)[None, :] self._moments = xe
