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.
python类inv()的实例源码
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 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 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)
test_linalg.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 25
收藏 0
点赞 0
评论 0
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))
test_defmatrix.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 36
收藏 0
点赞 0
评论 0
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)))
test_defmatrix.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 34
收藏 0
点赞 0
评论 0
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 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 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 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 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 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