mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-06-02 14:45:15 +02:00
remo0ved slices from models
slices are now handles by special indexing kern parts, such as coregionalisation, independent_outputs. The old slicing functionality has been removed simply to clean up the code a little. Now that input_slices still exist (and will continue to be useful) in kern.py. They do need a little work though, for the psi-statistics
This commit is contained in:
James Hensman
parent
ac842d51e6
commit
52ba8e4ba3
7 changed files with 103 additions and 175 deletions
152
GPy/kern/kern.py
152
GPy/kern/kern.py
|
|
@ -13,15 +13,9 @@ from prod import prod
|
|||
class kern(parameterised):
|
||||
def __init__(self, D, parts=[], input_slices=None):
|
||||
"""
|
||||
This kernel does 'compound' structures.
|
||||
This is the main kernel class for GPy. It handles multiple (additive) kernel functions, and keeps track of variaous things like which parameters live where.
|
||||
|
||||
The compund structure enables many features of GPy, including
|
||||
- Hierarchical models
|
||||
- Correleated output models
|
||||
- multi-view learning
|
||||
|
||||
Hadamard product and outer-product kernels will require a new class.
|
||||
This feature is currently WONTFIX. for small number sof inputs, you can use the sympy kernel for this.
|
||||
The technical code for kernels is divided into _parts_ (see e.g. rbf.py). This obnject contains a list of parts, which are computed additively. For multiplication, special _prod_ parts are used.
|
||||
|
||||
:param D: The dimensioality of the kernel's input space
|
||||
:type D: int
|
||||
|
|
@ -94,34 +88,6 @@ class kern(parameterised):
|
|||
self.param_slices.append(slice(count, count + p.Nparam))
|
||||
count += p.Nparam
|
||||
|
||||
def _process_slices(self, slices1=None, slices2=None):
|
||||
"""
|
||||
Format the slices so that they can easily be used.
|
||||
Both slices can be any of three things:
|
||||
- If None, the new points covary through every kernel part (default)
|
||||
- If a list of slices, the i^th slice specifies which data are affected by the i^th kernel part
|
||||
- If a list of booleans, specifying which kernel parts are active
|
||||
|
||||
if the second arg is False, return only slices1
|
||||
|
||||
returns actual lists of slice objects
|
||||
"""
|
||||
if slices1 is None:
|
||||
slices1 = [slice(None)] * self.Nparts
|
||||
elif all([type(s_i) is bool for s_i in slices1]):
|
||||
slices1 = [slice(None) if s_i else slice(0) for s_i in slices1]
|
||||
else:
|
||||
assert all([type(s_i) is slice for s_i in slices1]), "invalid slice objects"
|
||||
if slices2 is None:
|
||||
slices2 = [slice(None)] * self.Nparts
|
||||
elif slices2 is False:
|
||||
return slices1
|
||||
elif all([type(s_i) is bool for s_i in slices2]):
|
||||
slices2 = [slice(None) if s_i else slice(0) for s_i in slices2]
|
||||
else:
|
||||
assert all([type(s_i) is slice for s_i in slices2]), "invalid slice objects"
|
||||
return slices1, slices2
|
||||
|
||||
def __add__(self, other):
|
||||
assert self.D == other.D
|
||||
newkern = kern(self.D, self.parts + other.parts, self.input_slices + other.input_slices)
|
||||
|
|
@ -142,7 +108,7 @@ class kern(parameterised):
|
|||
:param other: the other kernel to be added
|
||||
:type other: GPy.kern
|
||||
"""
|
||||
return self +other
|
||||
return self + other
|
||||
|
||||
def add_orthogonal(self, other):
|
||||
"""
|
||||
|
|
@ -285,18 +251,19 @@ class kern(parameterised):
|
|||
|
||||
return sum([[name + '_' + n for n in k._get_param_names()] for name, k in zip(names, self.parts)], [])
|
||||
|
||||
def K(self, X, X2=None, slices1=None, slices2=None):
|
||||
def K(self, X, X2=None, which_parts='all'):
|
||||
if which_parts=='all':
|
||||
which_parts = [True]*self.Nparts
|
||||
assert X.shape[1] == self.D
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
if X2 is None:
|
||||
target = np.zeros((X.shape[0], X.shape[0]))
|
||||
[p.K(X[s1, i_s], None, target=target[s1, s2]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
[p.K(X[:, i_s], None, target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used]
|
||||
else:
|
||||
target = np.zeros((X.shape[0], X2.shape[0]))
|
||||
[p.K(X[s1, i_s], X2[s2, i_s], target=target[s1, s2]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
[p.K(X[:, i_s], X2[:,i_s], target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used]
|
||||
return target
|
||||
|
||||
def dK_dtheta(self, dL_dK, X, X2=None, slices1=None, slices2=None):
|
||||
def dK_dtheta(self, dL_dK, X, X2=None):
|
||||
"""
|
||||
:param dL_dK: An array of dL_dK derivaties, dL_dK
|
||||
:type dL_dK: Np.ndarray (N x M)
|
||||
|
|
@ -304,109 +271,94 @@ class kern(parameterised):
|
|||
:type X: np.ndarray (N x D)
|
||||
:param X2: Observed dara inputs (optional, defaults to X)
|
||||
:type X2: np.ndarray (M x D)
|
||||
:param slices1: a slice object for each kernel part, describing which data are affected by each kernel part
|
||||
:type slices1: list of slice objects, or list of booleans
|
||||
:param slices2: slices for X2
|
||||
"""
|
||||
assert X.shape[1] == self.D
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
target = np.zeros(self.Nparam)
|
||||
if X2 is None:
|
||||
[p.dK_dtheta(dL_dK[s1, s2], X[s1, i_s], None, target[ps]) for p, i_s, ps, s1, s2 in zip(self.parts, self.input_slices,self.param_slices, slices1, slices2)]
|
||||
[p.dK_dtheta(dL_dK, X[:, i_s], None, target[ps]) for p, i_s, ps, in zip(self.parts, self.input_slices, self.param_slices)]
|
||||
else:
|
||||
[p.dK_dtheta(dL_dK[s1, s2], X[s1, i_s], X2[s2, i_s], target[ps]) for p, i_s, ps, s1, s2 in zip(self.parts, self.input_slices,self.param_slices, slices1, slices2)]
|
||||
|
||||
[p.dK_dtheta(dL_dK, X[:, i_s], X2[:, i_s], target[ps]) for p, i_s, ps, in zip(self.parts, self.input_slices, self.param_slices)]
|
||||
|
||||
return self._transform_gradients(target)
|
||||
|
||||
def dK_dX(self, dL_dK, X, X2=None, slices1=None, slices2=None):
|
||||
def dK_dX(self, dL_dK, X, X2=None):
|
||||
if X2 is None:
|
||||
X2 = X
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
target = np.zeros_like(X)
|
||||
[p.dK_dX(dL_dK[s1, s2], X[s1, i_s], X2[s2, i_s], target[s1, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
if X2 is None:
|
||||
[p.dK_dX(dL_dK, X[:, i_s], None, target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
else:
|
||||
[p.dK_dX(dL_dK, X[:, i_s], X2[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target
|
||||
|
||||
def Kdiag(self, X, slices=None):
|
||||
def Kdiag(self, X, which_parts='all'):
|
||||
if which_parts=='all':
|
||||
which_parts = [True]*self.Nparts
|
||||
assert X.shape[1] == self.D
|
||||
slices = self._process_slices(slices, False)
|
||||
target = np.zeros(X.shape[0])
|
||||
[p.Kdiag(X[s, i_s], target=target[s]) for p, i_s, s in zip(self.parts, self.input_slices, slices)]
|
||||
[p.Kdiag(X[:, i_s], target=target) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target
|
||||
|
||||
def dKdiag_dtheta(self, dL_dKdiag, X, slices=None):
|
||||
def dKdiag_dtheta(self, dL_dKdiag, X):
|
||||
assert X.shape[1] == self.D
|
||||
assert len(dL_dKdiag.shape) == 1
|
||||
assert dL_dKdiag.size == X.shape[0]
|
||||
slices = self._process_slices(slices, False)
|
||||
target = np.zeros(self.Nparam)
|
||||
[p.dKdiag_dtheta(dL_dKdiag[s], X[s, i_s], target[ps]) for p, i_s, s, ps in zip(self.parts, self.input_slices, slices, self.param_slices)]
|
||||
[p.dKdiag_dtheta(dL_dKdiag, X[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)]
|
||||
return self._transform_gradients(target)
|
||||
|
||||
def dKdiag_dX(self, dL_dKdiag, X, slices=None):
|
||||
def dKdiag_dX(self, dL_dKdiag, X):
|
||||
assert X.shape[1] == self.D
|
||||
slices = self._process_slices(slices, False)
|
||||
target = np.zeros_like(X)
|
||||
[p.dKdiag_dX(dL_dKdiag[s], X[s, i_s], target[s, i_s]) for p, i_s, s in zip(self.parts, self.input_slices, slices)]
|
||||
[p.dKdiag_dX(dL_dKdiag, X[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target
|
||||
|
||||
def psi0(self, Z, mu, S, slices=None):
|
||||
slices = self._process_slices(slices, False)
|
||||
def psi0(self, Z, mu, S):
|
||||
target = np.zeros(mu.shape[0])
|
||||
[p.psi0(Z, mu[s], S[s], target[s]) for p, s in zip(self.parts, slices)]
|
||||
[p.psi0(Z[:,i_s], mu[:,i_s], S[:,i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target
|
||||
|
||||
def dpsi0_dtheta(self, dL_dpsi0, Z, mu, S, slices=None):
|
||||
slices = self._process_slices(slices, False)
|
||||
def dpsi0_dtheta(self, dL_dpsi0, Z, mu, S):
|
||||
target = np.zeros(self.Nparam)
|
||||
[p.dpsi0_dtheta(dL_dpsi0[s], Z, mu[s], S[s], target[ps]) for p, ps, s in zip(self.parts, self.param_slices, slices)]
|
||||
[p.dpsi0_dtheta(dL_dpsi0, Z[:,i_s], mu[:,i_s], S[:,i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)]
|
||||
return self._transform_gradients(target)
|
||||
|
||||
def dpsi0_dmuS(self, dL_dpsi0, Z, mu, S, slices=None):
|
||||
slices = self._process_slices(slices, False)
|
||||
def dpsi0_dmuS(self, dL_dpsi0, Z, mu, S):
|
||||
target_mu, target_S = np.zeros_like(mu), np.zeros_like(S)
|
||||
[p.dpsi0_dmuS(dL_dpsi0, Z, mu[s], S[s], target_mu[s], target_S[s]) for p, s in zip(self.parts, slices)]
|
||||
[p.dpsi0_dmuS(dL_dpsi0, Z[:,i_s], mu[:,i_s], S[:,i_s], target_mu[:,i_s], target_S[:,i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target_mu, target_S
|
||||
|
||||
def psi1(self, Z, mu, S, slices1=None, slices2=None):
|
||||
"""Think N,M,Q """
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
def psi1(self, Z, mu, S):
|
||||
target = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
[p.psi1(Z[s2], mu[s1], S[s1], target[s1, s2]) for p, s1, s2 in zip(self.parts, slices1, slices2)]
|
||||
[p.psi1(Z[:,i_s], mu[:,i_s], S[:,i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target
|
||||
|
||||
def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S, slices1=None, slices2=None):
|
||||
"""N,M,(Ntheta)"""
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S):
|
||||
target = np.zeros((self.Nparam))
|
||||
[p.dpsi1_dtheta(dL_dpsi1[s2, s1], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[ps]) for p, ps, s1, s2, i_s in zip(self.parts, self.param_slices, slices1, slices2, self.input_slices)]
|
||||
[p.dpsi1_dtheta(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)]
|
||||
return self._transform_gradients(target)
|
||||
|
||||
def dpsi1_dZ(self, dL_dpsi1, Z, mu, S, slices1=None, slices2=None):
|
||||
"""N,M,Q"""
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
def dpsi1_dZ(self, dL_dpsi1, Z, mu, S):
|
||||
target = np.zeros_like(Z)
|
||||
[p.dpsi1_dZ(dL_dpsi1[s2, s1], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[s2, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
[p.dpsi1_dZ(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target
|
||||
|
||||
def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S, slices1=None, slices2=None):
|
||||
def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
|
||||
"""return shapes are N,M,Q"""
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
|
||||
[p.dpsi1_dmuS(dL_dpsi1[s2, s1], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target_mu[s1, i_s], target_S[s1, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
[p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target_mu, target_S
|
||||
|
||||
def psi2(self, Z, mu, S, slices1=None, slices2=None):
|
||||
def psi2(self, Z, mu, S):
|
||||
"""
|
||||
:param Z: np.ndarray of inducing inputs (M x Q)
|
||||
:param mu, S: np.ndarrays of means and variances (each N x Q)
|
||||
:returns psi2: np.ndarray (N,M,M)
|
||||
"""
|
||||
target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
[p.psi2(Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[s1, s2, s2]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
[p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
|
||||
# compute the "cross" terms
|
||||
#TODO: input_slices needed
|
||||
for p1, p2 in itertools.combinations(self.parts, 2):
|
||||
# white doesn;t combine with anything
|
||||
if p1.name == 'white' or p2.name == 'white':
|
||||
|
|
@ -434,14 +386,12 @@ class kern(parameterised):
|
|||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
return target
|
||||
|
||||
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, slices1=None, slices2=None):
|
||||
"""Returns shape (N,M,M,Ntheta)"""
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S):
|
||||
target = np.zeros(self.Nparam)
|
||||
[p.dpsi2_dtheta(dL_dpsi2[s1, s2, s2], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[ps]) for p, i_s, s1, s2, ps in zip(self.parts, self.input_slices, slices1, slices2, self.param_slices)]
|
||||
[p.dpsi2_dtheta(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)]
|
||||
|
||||
# compute the "cross" terms
|
||||
# TODO: better looping
|
||||
# TODO: better looping, input_slices
|
||||
for i1, i2 in itertools.combinations(range(len(self.parts)), 2):
|
||||
p1, p2 = self.parts[i1], self.parts[i2]
|
||||
# ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2]
|
||||
|
|
@ -478,12 +428,12 @@ class kern(parameterised):
|
|||
|
||||
return self._transform_gradients(target)
|
||||
|
||||
def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, slices1=None, slices2=None):
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
def dpsi2_dZ(self, dL_dpsi2, Z, mu, S):
|
||||
target = np.zeros_like(Z)
|
||||
[p.dpsi2_dZ(dL_dpsi2[s1, s2, s2], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[s2, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
[p.dpsi2_dZ(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
|
||||
# compute the "cross" terms
|
||||
#TODO: we need input_slices here.
|
||||
for p1, p2 in itertools.combinations(self.parts, 2):
|
||||
# white doesn;t combine with anything
|
||||
if p1.name == 'white' or p2.name == 'white':
|
||||
|
|
@ -506,16 +456,14 @@ class kern(parameterised):
|
|||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
|
||||
|
||||
return target * 2.
|
||||
|
||||
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, slices1=None, slices2=None):
|
||||
"""return shapes are N,M,M,Q"""
|
||||
slices1, slices2 = self._process_slices(slices1, slices2)
|
||||
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S):
|
||||
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
|
||||
[p.dpsi2_dmuS(dL_dpsi2[s1, s2, s2], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target_mu[s1, i_s], target_S[s1, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)]
|
||||
[p.dpsi2_dmuS(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
|
||||
# compute the "cross" terms
|
||||
#TODO: we need input_slices here.
|
||||
for p1, p2 in itertools.combinations(self.parts, 2):
|
||||
# white doesn;t combine with anything
|
||||
if p1.name == 'white' or p2.name == 'white':
|
||||
|
|
|
|||
|
|
@ -19,7 +19,6 @@ class GP(model):
|
|||
:parm likelihood: a GPy likelihood
|
||||
:param normalize_X: whether to normalize the input data before computing (predictions will be in original scales)
|
||||
:type normalize_X: False|True
|
||||
:param Xslices: how the X,Y data co-vary in the kernel (i.e. which "outputs" they correspond to). See (link:slicing)
|
||||
:rtype: model object
|
||||
:param epsilon_ep: convergence criterion for the Expectation Propagation algorithm, defaults to 0.1
|
||||
:param powerep: power-EP parameters [$\eta$,$\delta$], defaults to [1.,1.]
|
||||
|
|
@ -28,10 +27,9 @@ class GP(model):
|
|||
.. Note:: Multiple independent outputs are allowed using columns of Y
|
||||
|
||||
"""
|
||||
def __init__(self, X, likelihood, kernel, normalize_X=False, Xslices=None):
|
||||
def __init__(self, X, likelihood, kernel, normalize_X=False):
|
||||
|
||||
# parse arguments
|
||||
self.Xslices = Xslices
|
||||
self.X = X
|
||||
assert len(self.X.shape) == 2
|
||||
self.N, self.Q = self.X.shape
|
||||
|
|
@ -64,12 +62,12 @@ class GP(model):
|
|||
return np.zeros_like(self.Z)
|
||||
|
||||
def _set_params(self, p):
|
||||
self.kern._set_params_transformed(p[:self.kern.Nparam])
|
||||
self.kern._set_params_transformed(p[:self.kern.Nparam_transformed()])
|
||||
# self.likelihood._set_params(p[self.kern.Nparam:]) # test by Nicolas
|
||||
self.likelihood._set_params(p[self.kern.Nparam_transformed():]) # test by Nicolas
|
||||
|
||||
|
||||
self.K = self.kern.K(self.X, slices1=self.Xslices, slices2=self.Xslices)
|
||||
self.K = self.kern.K(self.X)
|
||||
self.K += self.likelihood.covariance_matrix
|
||||
|
||||
self.Ki, self.L, self.Li, self.K_logdet = pdinv(self.K)
|
||||
|
|
@ -92,7 +90,7 @@ class GP(model):
|
|||
"""
|
||||
Approximates a non-gaussian likelihood using Expectation Propagation
|
||||
|
||||
For a Gaussian (or direct: TODO) likelihood, no iteration is required:
|
||||
For a Gaussian likelihood, no iteration is required:
|
||||
this function does nothing
|
||||
"""
|
||||
self.likelihood.fit_full(self.kern.K(self.X))
|
||||
|
|
@ -122,31 +120,33 @@ class GP(model):
|
|||
"""
|
||||
The gradient of all parameters.
|
||||
|
||||
For the kernel parameters, use the chain rule via dL_dK
|
||||
|
||||
For the likelihood parameters, pass in alpha = K^-1 y
|
||||
Note, we use the chain rule: dL_dtheta = dL_dK * d_K_dtheta
|
||||
"""
|
||||
return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X, slices1=self.Xslices, slices2=self.Xslices), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
|
||||
return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
|
||||
|
||||
def _raw_predict(self, _Xnew, slices=None, full_cov=False):
|
||||
def _raw_predict(self, _Xnew, which_parts='all', full_cov=False):
|
||||
"""
|
||||
Internal helper function for making predictions, does not account
|
||||
for normalization or likelihood
|
||||
|
||||
#TODO: which_parts does nothing
|
||||
|
||||
|
||||
"""
|
||||
Kx = self.kern.K(self.X, _Xnew, slices1=self.Xslices, slices2=slices)
|
||||
Kx = self.kern.K(self.X, _Xnew,which_parts=which_parts)
|
||||
mu = np.dot(np.dot(Kx.T, self.Ki), self.likelihood.Y)
|
||||
KiKx = np.dot(self.Ki, Kx)
|
||||
if full_cov:
|
||||
Kxx = self.kern.K(_Xnew, slices1=slices, slices2=slices)
|
||||
Kxx = self.kern.K(_Xnew, which_parts=which_parts)
|
||||
var = Kxx - np.dot(KiKx.T, Kx)
|
||||
else:
|
||||
Kxx = self.kern.Kdiag(_Xnew, slices=slices)
|
||||
Kxx = self.kern.Kdiag(_Xnew, which_parts=which_parts)
|
||||
var = Kxx - np.sum(np.multiply(KiKx, Kx), 0)
|
||||
var = var[:, None]
|
||||
return mu, var
|
||||
|
||||
|
||||
def predict(self, Xnew, slices=None, full_cov=False):
|
||||
def predict(self, Xnew, which_parts='all', full_cov=False):
|
||||
"""
|
||||
Predict the function(s) at the new point(s) Xnew.
|
||||
|
||||
|
|
@ -154,19 +154,14 @@ class GP(model):
|
|||
---------
|
||||
:param Xnew: The points at which to make a prediction
|
||||
:type Xnew: np.ndarray, Nnew x self.Q
|
||||
:param slices: specifies which outputs kernel(s) the Xnew correspond to (see below)
|
||||
:type slices: (None, list of slice objects, list of ints)
|
||||
:param which_parts: specifies which outputs kernel(s) to use in prediction
|
||||
:type which_parts: ('all', list of bools)
|
||||
:param full_cov: whether to return the folll covariance matrix, or just the diagonal
|
||||
:type full_cov: bool
|
||||
:rtype: posterior mean, a Numpy array, Nnew x self.D
|
||||
:rtype: posterior variance, a Numpy array, Nnew x 1 if full_cov=False, Nnew x Nnew otherwise
|
||||
:rtype: lower and upper boundaries of the 95% confidence intervals, Numpy arrays, Nnew x self.D
|
||||
|
||||
.. Note:: "slices" specifies how the the points X_new co-vary wich the training points.
|
||||
|
||||
- If None, the new points covary throigh every kernel part (default)
|
||||
- If a list of slices, the i^th slice specifies which data are affected by the i^th kernel part
|
||||
- If a list of booleans, specifying which kernel parts are active
|
||||
|
||||
If full_cov and self.D > 1, the return shape of var is Nnew x Nnew x self.D. If self.D == 1, the return shape is Nnew x Nnew.
|
||||
This is to allow for different normalizations of the output dimensions.
|
||||
|
|
@ -174,15 +169,15 @@ class GP(model):
|
|||
"""
|
||||
# normalize X values
|
||||
Xnew = (Xnew.copy() - self._Xmean) / self._Xstd
|
||||
mu, var = self._raw_predict(Xnew, slices, full_cov)
|
||||
mu, var = self._raw_predict(Xnew, which_parts, full_cov)
|
||||
|
||||
# now push through likelihood TODO
|
||||
# now push through likelihood
|
||||
mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov)
|
||||
|
||||
return mean, var, _025pm, _975pm
|
||||
|
||||
|
||||
def plot_f(self, samples=0, plot_limits=None, which_data='all', which_functions='all', resolution=None, full_cov=False):
|
||||
def plot_f(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False):
|
||||
"""
|
||||
Plot the GP's view of the world, where the data is normalized and the likelihood is Gaussian
|
||||
|
||||
|
|
@ -190,8 +185,8 @@ class GP(model):
|
|||
:param which_data: which if the training data to plot (default all)
|
||||
:type which_data: 'all' or a slice object to slice self.X, self.Y
|
||||
:param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
|
||||
:param which_functions: which of the kernel functions to plot (additively)
|
||||
:type which_functions: list of bools
|
||||
:param which_parts: which of the kernel functions to plot (additively)
|
||||
:type which_parts: 'all', or list of bools
|
||||
:param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
|
||||
|
||||
Plot the posterior of the GP.
|
||||
|
|
@ -202,19 +197,17 @@ class GP(model):
|
|||
Can plot only part of the data and part of the posterior functions using which_data and which_functions
|
||||
Plot the data's view of the world, with non-normalized values and GP predictions passed through the likelihood
|
||||
"""
|
||||
if which_functions == 'all':
|
||||
which_functions = [True] * self.kern.Nparts
|
||||
if which_data == 'all':
|
||||
which_data = slice(None)
|
||||
|
||||
if self.X.shape[1] == 1:
|
||||
Xnew, xmin, xmax = x_frame1D(self.X, plot_limits=plot_limits)
|
||||
if samples == 0:
|
||||
m, v = self._raw_predict(Xnew, slices=which_functions)
|
||||
m, v = self._raw_predict(Xnew, which_parts=which_parts)
|
||||
gpplot(Xnew, m, m - 2 * np.sqrt(v), m + 2 * np.sqrt(v))
|
||||
pb.plot(self.X[which_data], self.likelihood.Y[which_data], 'kx', mew=1.5)
|
||||
else:
|
||||
m, v = self._raw_predict(Xnew, slices=which_functions, full_cov=True)
|
||||
m, v = self._raw_predict(Xnew, which_parts=which_parts, full_cov=True)
|
||||
Ysim = np.random.multivariate_normal(m.flatten(), v, samples)
|
||||
gpplot(Xnew, m, m - 2 * np.sqrt(np.diag(v)[:, None]), m + 2 * np.sqrt(np.diag(v))[:, None])
|
||||
for i in range(samples):
|
||||
|
|
@ -230,7 +223,7 @@ class GP(model):
|
|||
elif self.X.shape[1] == 2:
|
||||
resolution = resolution or 50
|
||||
Xnew, xmin, xmax, xx, yy = x_frame2D(self.X, plot_limits, resolution)
|
||||
m, v = self._raw_predict(Xnew, slices=which_functions)
|
||||
m, v = self._raw_predict(Xnew, which_parts=which_parts)
|
||||
m = m.reshape(resolution, resolution).T
|
||||
pb.contour(xx, yy, m, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
|
||||
pb.scatter(Xorig[:, 0], Xorig[:, 1], 40, Yorig, linewidth=0, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max())
|
||||
|
|
@ -246,8 +239,6 @@ class GP(model):
|
|||
|
||||
"""
|
||||
# TODO include samples
|
||||
if which_functions == 'all':
|
||||
which_functions = [True] * self.kern.Nparts
|
||||
if which_data == 'all':
|
||||
which_data = slice(None)
|
||||
|
||||
|
|
@ -256,7 +247,7 @@ class GP(model):
|
|||
Xu = self.X * self._Xstd + self._Xmean # NOTE self.X are the normalized values now
|
||||
|
||||
Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits)
|
||||
m, var, lower, upper = self.predict(Xnew, slices=which_functions)
|
||||
m, var, lower, upper = self.predict(Xnew, which_parts=which_parts)
|
||||
gpplot(Xnew, m, lower, upper)
|
||||
pb.plot(Xu[which_data], self.likelihood.data[which_data], 'kx', mew=1.5)
|
||||
if self.has_uncertain_inputs:
|
||||
|
|
@ -277,7 +268,7 @@ class GP(model):
|
|||
resolution = resolution or 50
|
||||
Xnew, xx, yy, xmin, xmax = x_frame2D(self.X, plot_limits, resolution)
|
||||
x, y = np.linspace(xmin[0], xmax[0], resolution), np.linspace(xmin[1], xmax[1], resolution)
|
||||
m, var, lower, upper = self.predict(Xnew, slices=which_functions)
|
||||
m, var, lower, upper = self.predict(Xnew, which_parts=which_parts)
|
||||
m = m.reshape(resolution, resolution).T
|
||||
pb.contour(x, y, m, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
|
||||
Yf = self.likelihood.Y.flatten()
|
||||
|
|
|
|||
|
|
@ -11,26 +11,24 @@ class GP_regression(GP):
|
|||
"""
|
||||
Gaussian Process model for regression
|
||||
|
||||
This is a thin wrapper around the GP class, with a set of sensible defalts
|
||||
This is a thin wrapper around the models.GP class, with a set of sensible defalts
|
||||
|
||||
:param X: input observations
|
||||
:param Y: observed values
|
||||
:param kernel: a GPy kernel, defaults to rbf+white
|
||||
:param kernel: a GPy kernel, defaults to rbf
|
||||
:param normalize_X: whether to normalize the input data before computing (predictions will be in original scales)
|
||||
:type normalize_X: False|True
|
||||
:param normalize_Y: whether to normalize the input data before computing (predictions will be in original scales)
|
||||
:type normalize_Y: False|True
|
||||
:param Xslices: how the X,Y data co-vary in the kernel (i.e. which "outputs" they correspond to). See (link:slicing)
|
||||
:rtype: model object
|
||||
|
||||
.. Note:: Multiple independent outputs are allowed using columns of Y
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self,X,Y,kernel=None,normalize_X=False,normalize_Y=False, Xslices=None):
|
||||
def __init__(self,X,Y,kernel=None,normalize_X=False,normalize_Y=False):
|
||||
if kernel is None:
|
||||
kernel = kern.rbf(X.shape[1])
|
||||
|
||||
likelihood = likelihoods.Gaussian(Y,normalize=normalize_Y)
|
||||
|
||||
GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X, Xslices=Xslices)
|
||||
GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X)
|
||||
|
|
|
|||
|
|
@ -23,20 +23,19 @@ class generalized_FITC(sparse_GP):
|
|||
:type X_variance: np.ndarray (N x Q) | None
|
||||
:param Z: inducing inputs (optional, see note)
|
||||
:type Z: np.ndarray (M x Q) | None
|
||||
:param Zslices: slices for the inducing inputs (see slicing TODO: link)
|
||||
:param M : Number of inducing points (optional, default 10. Ignored if Z is not None)
|
||||
:type M: int
|
||||
:param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
|
||||
:type normalize_(X|Y): bool
|
||||
"""
|
||||
|
||||
def __init__(self, X, likelihood, kernel, Z, X_variance=None, Xslices=None,Zslices=None, normalize_X=False):
|
||||
def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False):
|
||||
|
||||
self.Z = Z
|
||||
self.M = self.Z.shape[0]
|
||||
self._precision = likelihood.precision
|
||||
|
||||
sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_variance=None, Xslices=None,Zslices=None, normalize_X=False)
|
||||
sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_variance=None, normalize_X=False)
|
||||
|
||||
def _set_params(self, p):
|
||||
self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
|
||||
|
|
@ -145,7 +144,7 @@ class generalized_FITC(sparse_GP):
|
|||
D = 0.5*np.trace(self.Cpsi1VVpsi1)
|
||||
return A+C+D
|
||||
|
||||
def _raw_predict(self, Xnew, slices, full_cov=False):
|
||||
def _raw_predict(self, Xnew, which_parts, full_cov=False):
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
"""
|
||||
Make a prediction for the generalized FITC model
|
||||
|
|
@ -174,16 +173,16 @@ class generalized_FITC(sparse_GP):
|
|||
self.mu_H = mu_H
|
||||
Sigma_H = C + np.dot(mu_u,np.dot(self.Sigma,mu_u.T))
|
||||
# q(f_star|y) = N(f_star|mu_star,sigma2_star)
|
||||
Kx = self.kern.K(self.Z, Xnew)
|
||||
Kx = self.kern.K(self.Z, Xnew, which_parts=which_parts)
|
||||
KR0T = np.dot(Kx.T,self.Lmi.T)
|
||||
mu_star = np.dot(KR0T,mu_H)
|
||||
if full_cov:
|
||||
Kxx = self.kern.K(Xnew)
|
||||
Kxx = self.kern.K(Xnew,which_parts=which_parts)
|
||||
var = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T))
|
||||
else:
|
||||
Kxx = self.kern.Kdiag(Xnew)
|
||||
Kxx_ = self.kern.K(Xnew)
|
||||
var_ = Kxx_ + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T))
|
||||
Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts)
|
||||
Kxx_ = self.kern.K(Xnew,which_parts=which_parts) # TODO: RA, is this line needed?
|
||||
var_ = Kxx_ + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) # TODO: RA, is this line needed?
|
||||
var = (Kxx + np.sum(KR0T.T*np.dot(Sigma_H - np.eye(self.M),KR0T.T),0))[:,None]
|
||||
return mu_star[:,None],var
|
||||
else:
|
||||
|
|
|
|||
|
|
@ -9,10 +9,6 @@ from .. import kern
|
|||
from GP import GP
|
||||
from scipy import linalg
|
||||
|
||||
#Still TODO:
|
||||
# make use of slices properly (kernel can now do this)
|
||||
# enable heteroscedatic noise (kernel will need to compute psi2 as a (NxMxM) array)
|
||||
|
||||
class sparse_GP(GP):
|
||||
"""
|
||||
Variational sparse GP model
|
||||
|
|
@ -27,19 +23,16 @@ class sparse_GP(GP):
|
|||
:type X_variance: np.ndarray (N x Q) | None
|
||||
:param Z: inducing inputs (optional, see note)
|
||||
:type Z: np.ndarray (M x Q) | None
|
||||
:param Zslices: slices for the inducing inputs (see slicing TODO: link)
|
||||
:param M : Number of inducing points (optional, default 10. Ignored if Z is not None)
|
||||
:type M: int
|
||||
:param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
|
||||
:type normalize_(X|Y): bool
|
||||
"""
|
||||
|
||||
def __init__(self, X, likelihood, kernel, Z, X_variance=None, Xslices=None,Zslices=None, normalize_X=False):
|
||||
def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False):
|
||||
self.scale_factor = 100.0# a scaling factor to help keep the algorithm stable
|
||||
self.auto_scale_factor = False
|
||||
self.Z = Z
|
||||
self.Zslices = Zslices
|
||||
self.Xslices = Xslices
|
||||
self.M = Z.shape[0]
|
||||
self.likelihood = likelihood
|
||||
|
||||
|
|
@ -50,7 +43,7 @@ class sparse_GP(GP):
|
|||
self.has_uncertain_inputs=True
|
||||
self.X_variance = X_variance
|
||||
|
||||
GP.__init__(self, X, likelihood, kernel=kernel, normalize_X=normalize_X, Xslices=Xslices)
|
||||
GP.__init__(self, X, likelihood, kernel=kernel, normalize_X=normalize_X)
|
||||
|
||||
#normalize X uncertainty also
|
||||
if self.has_uncertain_inputs:
|
||||
|
|
@ -65,13 +58,12 @@ class sparse_GP(GP):
|
|||
self.psi1 = self.kern.psi1(self.Z,self.X, self.X_variance).T
|
||||
self.psi2 = self.kern.psi2(self.Z,self.X, self.X_variance)
|
||||
else:
|
||||
self.psi0 = self.kern.Kdiag(self.X,slices=self.Xslices)
|
||||
self.psi0 = self.kern.Kdiag(self.X)
|
||||
self.psi1 = self.kern.K(self.Z,self.X)
|
||||
self.psi2 = None
|
||||
|
||||
def _computations(self):
|
||||
#TODO: find routine to multiply triangular matrices
|
||||
#TODO: slices for psi statistics (easy enough)
|
||||
|
||||
sf = self.scale_factor
|
||||
sf2 = sf**2
|
||||
|
|
@ -252,16 +244,16 @@ class sparse_GP(GP):
|
|||
dL_dZ += self.kern.dK_dX(self.dL_dpsi1,self.Z,self.X)
|
||||
return dL_dZ
|
||||
|
||||
def _raw_predict(self, Xnew, slices, full_cov=False):
|
||||
def _raw_predict(self, Xnew, which_parts='all', full_cov=False):
|
||||
"""Internal helper function for making predictions, does not account for normalization"""
|
||||
|
||||
Kx = self.kern.K(self.Z, Xnew)
|
||||
mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V)
|
||||
if full_cov:
|
||||
Kxx = self.kern.K(Xnew)
|
||||
Kxx = self.kern.K(Xnew,which_parts=which_parts)
|
||||
var = Kxx - mdot(Kx.T, (self.Kmmi - self.C/self.scale_factor**2), Kx) #NOTE this won't work for plotting
|
||||
else:
|
||||
Kxx = self.kern.Kdiag(Xnew)
|
||||
Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts)
|
||||
var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0)
|
||||
|
||||
return mu,var[:,None]
|
||||
|
|
|
|||
|
|
@ -13,7 +13,7 @@ class sparse_GP_regression(sparse_GP):
|
|||
"""
|
||||
Gaussian Process model for regression
|
||||
|
||||
This is a thin wrapper around the GP class, with a set of sensible defalts
|
||||
This is a thin wrapper around the sparse_GP class, with a set of sensible defalts
|
||||
|
||||
:param X: input observations
|
||||
:param Y: observed values
|
||||
|
|
@ -22,25 +22,25 @@ class sparse_GP_regression(sparse_GP):
|
|||
:type normalize_X: False|True
|
||||
:param normalize_Y: whether to normalize the input data before computing (predictions will be in original scales)
|
||||
:type normalize_Y: False|True
|
||||
:param Xslices: how the X,Y data co-vary in the kernel (i.e. which "outputs" they correspond to). See (link:slicing)
|
||||
:rtype: model object
|
||||
|
||||
.. Note:: Multiple independent outputs are allowed using columns of Y
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self,X,Y,kernel=None,normalize_X=False,normalize_Y=False, Xslices=None,Z=None, M=10):
|
||||
#kern defaults to rbf
|
||||
def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False, Z=None, M=10):
|
||||
#kern defaults to rbf (plus white for stability)
|
||||
if kernel is None:
|
||||
kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1],1e-3)
|
||||
|
||||
#Z defaults to a subset of the data
|
||||
if Z is None:
|
||||
Z = np.random.permutation(X.copy())[:M]
|
||||
i = np.random.permutation(X.shape[0])[:M]
|
||||
Z = X[i].copy()
|
||||
else:
|
||||
assert Z.shape[1]==X.shape[1]
|
||||
|
||||
#likelihood defaults to Gaussian
|
||||
likelihood = likelihoods.Gaussian(Y,normalize=normalize_Y)
|
||||
|
||||
sparse_GP.__init__(self, X, likelihood, kernel, Z, normalize_X=normalize_X, Xslices=Xslices)
|
||||
sparse_GP.__init__(self, X, likelihood, kernel, Z, normalize_X=normalize_X)
|
||||
|
|
|
|||
|
|
@ -14,7 +14,7 @@ from .. import likelihoods
|
|||
from .. import kern
|
||||
|
||||
class warpedGP(GP):
|
||||
def __init__(self, X, Y, kernel=None, warping_function = None, warping_terms = 3, normalize_X=False, normalize_Y=False, Xslices=None):
|
||||
def __init__(self, X, Y, kernel=None, warping_function = None, warping_terms = 3, normalize_X=False, normalize_Y=False):
|
||||
|
||||
if kernel is None:
|
||||
kernel = kern.rbf(X.shape[1])
|
||||
|
|
@ -28,7 +28,7 @@ class warpedGP(GP):
|
|||
self.predict_in_warped_space = False
|
||||
likelihood = likelihoods.Gaussian(self.transform_data(), normalize=normalize_Y)
|
||||
|
||||
GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X, Xslices=Xslices)
|
||||
GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X)
|
||||
|
||||
def _set_params(self, x):
|
||||
self.warping_params = x[:self.warping_function.num_parameters]
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue