apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-09 20:12:38 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'params' of github.com:SheffieldML/GPy into params

Browse source 
Conflicts:
	GPy/core/parameterization/param.py
	GPy/core/parameterization/parameter_core.py
	GPy/core/parameterization/parameterized.py
This commit is contained in:
Max Zwiessele 2014-02-10 15:21:09 +00:00
commit 6a068775f5
10 changed files with 404 additions and 301 deletions
Download patch file Download diff file Expand all files Collapse all files

77
GPy/core/parameterization/param.py
View file

@ -4,19 +4,19 @@
import itertools import itertools
import numpy import numpy
from parameter_core import Constrainable, adjust_name_for_printing from parameter_core import Constrainable, adjust_name_for_printing
from array_core import ObservableArray, ParamList from array_core import ObservableArray
###### printing ###### printing
__constraints_name__ = "Constraint" __constraints_name__ = "Constraint"
__index_name__ = "Index" __index_name__ = "Index"
__tie_name__ = "Tied to" __tie_name__ = "Tied to"
__precision__ = numpy.get_printoptions()['precision'] # numpy printing precision used, sublassing numpy ndarray after all __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision used, sublassing numpy ndarray after all
__print_threshold__ = 5 __print_threshold__ = 5
###### ######
class Float(numpy.float64, Constrainable): class Float(numpy.float64, Constrainable):
def __init__(self, f, base): def __init__(self, f, base):
super(Float, self).__init__(f) super(Float,self).__init__(f)
self._base = base self._base = base
@ -50,7 +50,7 @@ class Param(ObservableArray, Constrainable):
WARNING: This overrides the functionality of x==y!!! WARNING: This overrides the functionality of x==y!!!
Use numpy.equal(x,y) for element-wise equality testing. Use numpy.equal(x,y) for element-wise equality testing.
""" """
__array_priority__ = 0 # Never give back Param __array_priority__ = 0 # Never give back Param
_fixes_ = None _fixes_ = None
def __new__(cls, name, input_array, *args, **kwargs): def __new__(cls, name, input_array, *args, **kwargs):
obj = numpy.atleast_1d(super(Param, cls).__new__(cls, input_array=input_array)) obj = numpy.atleast_1d(super(Param, cls).__new__(cls, input_array=input_array))
@ -185,7 +185,7 @@ class Param(ObservableArray, Constrainable):
Note: For now only one parameter can have ties, so all of a parameter Note: For now only one parameter can have ties, so all of a parameter
will be removed, when re-tieing! will be removed, when re-tieing!
""" """
# Note: this method will tie to the parameter which is the last in #Note: this method will tie to the parameter which is the last in
# the chain of ties. Thus, if you tie to a tied parameter, # the chain of ties. Thus, if you tie to a tied parameter,
# this tie will be created to the parameter the param is tied # this tie will be created to the parameter the param is tied
# to. # to.
@ -197,7 +197,7 @@ class Param(ObservableArray, Constrainable):
try: try:
if self._original_: if self._original_:
self[:] = param self[:] = param
else: # this happens when indexing created a copy of the array else: # this happens when indexing created a copy of the array
self._direct_parent_._get_original(self)[self._current_slice_] = param self._direct_parent_._get_original(self)[self._current_slice_] = param
except ValueError: except ValueError:
raise ValueError("Trying to tie {} with shape {} to {} with shape {}".format(self.name, self.shape, param.name, param.shape)) raise ValueError("Trying to tie {} with shape {} to {} with shape {}".format(self.name, self.shape, param.name, param.shape))
@ -320,7 +320,7 @@ class Param(ObservableArray, Constrainable):
continue continue
if isinstance(si, slice): if isinstance(si, slice):
a = si.indices(self._realshape_[i])[0] a = si.indices(self._realshape_[i])[0]
elif isinstance(si, (list, numpy.ndarray, tuple)): elif isinstance(si, (list,numpy.ndarray,tuple)):
a = si[0] a = si[0]
else: a = si else: a = si
if a < 0: if a < 0:
@ -475,20 +475,20 @@ class ParamConcatenation(object):
self.params.append(p) self.params.append(p)
self._param_sizes = [p.size for p in self.params] self._param_sizes = [p.size for p in self.params]
startstops = numpy.cumsum([0] + self._param_sizes) startstops = numpy.cumsum([0] + self._param_sizes)
self._param_slices_ = [slice(start, stop) for start, stop in zip(startstops, startstops[1:])] self._param_slices_ = [slice(start, stop) for start,stop in zip(startstops, startstops[1:])]
#=========================================================================== #===========================================================================
# Get/set items, enable broadcasting # Get/set items, enable broadcasting
#=========================================================================== #===========================================================================
def __getitem__(self, s): def __getitem__(self, s):
ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True;
params = [p._get_params()[ind[ps]] for p, ps in zip(self.params, self._param_slices_) if numpy.any(p._get_params()[ind[ps]])] params = [p._get_params()[ind[ps]] for p,ps in zip(self.params, self._param_slices_) if numpy.any(p._get_params()[ind[ps]])]
if len(params) == 1: return params[0] if len(params)==1: return params[0]
return ParamConcatenation(params) return ParamConcatenation(params)
def __setitem__(self, s, val, update=True): def __setitem__(self, s, val, update=True):
ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True;
vals = self._vals(); vals[s] = val; del val vals = self._vals(); vals[s] = val; del val
[numpy.place(p, ind[ps], vals[ps]) and p._notify_tied_parameters() [numpy.place(p, ind[ps], vals[ps]) and p._notify_tied_parameters()
for p, ps in zip(self.params, self._param_slices_)] for p, ps in zip(self.params, self._param_slices_)]
if update: if update:
self.params[0]._highest_parent_.parameters_changed() self.params[0]._highest_parent_.parameters_changed()
def _vals(self): def _vals(self):
@ -496,38 +496,55 @@ class ParamConcatenation(object):
#=========================================================================== #===========================================================================
# parameter operations: # parameter operations:
#=========================================================================== #===========================================================================
def update_all_params(self):
self.params[0]._highest_parent_.parameters_changed()
def constrain(self, constraint, warning=True): def constrain(self, constraint, warning=True):
[param.constrain(constraint) for param in self.params] [param.constrain(constraint, update=False) for param in self.params]
self.update_all_params()
constrain.__doc__ = Param.constrain.__doc__ constrain.__doc__ = Param.constrain.__doc__
def constrain_positive(self, warning=True): def constrain_positive(self, warning=True):
[param.constrain_positive(warning) for param in self.params] [param.constrain_positive(warning, update=False) for param in self.params]
self.update_all_params()
constrain_positive.__doc__ = Param.constrain_positive.__doc__ constrain_positive.__doc__ = Param.constrain_positive.__doc__
def constrain_fixed(self, warning=True): def constrain_fixed(self, warning=True):
[param.constrain_fixed(warning) for param in self.params] [param.constrain_fixed(warning) for param in self.params]
constrain_fixed.__doc__ = Param.constrain_fixed.__doc__ constrain_fixed.__doc__ = Param.constrain_fixed.__doc__
fix = constrain_fixed fix = constrain_fixed
def constrain_negative(self, warning=True): def constrain_negative(self, warning=True):
[param.constrain_negative(warning) for param in self.params] [param.constrain_negative(warning, update=False) for param in self.params]
self.update_all_params()
constrain_negative.__doc__ = Param.constrain_negative.__doc__ constrain_negative.__doc__ = Param.constrain_negative.__doc__
def constrain_bounded(self, lower, upper, warning=True): def constrain_bounded(self, lower, upper, warning=True):
[param.constrain_bounded(lower, upper, warning) for param in self.params] [param.constrain_bounded(lower, upper, warning, update=False) for param in self.params]
self.update_all_params()
constrain_bounded.__doc__ = Param.constrain_bounded.__doc__ constrain_bounded.__doc__ = Param.constrain_bounded.__doc__
def unconstrain(self, *constraints): def unconstrain(self, *constraints):
[param.unconstrain(*constraints) for param in self.params] [param.unconstrain(*constraints) for param in self.params]
unconstrain.__doc__ = Param.unconstrain.__doc__ unconstrain.__doc__ = Param.unconstrain.__doc__
def unconstrain_negative(self): def unconstrain_negative(self):
[param.unconstrain_negative() for param in self.params] [param.unconstrain_negative() for param in self.params]
unconstrain_negative.__doc__ = Param.unconstrain_negative.__doc__ unconstrain_negative.__doc__ = Param.unconstrain_negative.__doc__
def unconstrain_positive(self): def unconstrain_positive(self):
[param.unconstrain_positive() for param in self.params] [param.unconstrain_positive() for param in self.params]
unconstrain_positive.__doc__ = Param.unconstrain_positive.__doc__ unconstrain_positive.__doc__ = Param.unconstrain_positive.__doc__
def unconstrain_fixed(self): def unconstrain_fixed(self):
[param.unconstrain_fixed() for param in self.params] [param.unconstrain_fixed() for param in self.params]
unconstrain_fixed.__doc__ = Param.unconstrain_fixed.__doc__ unconstrain_fixed.__doc__ = Param.unconstrain_fixed.__doc__
unfix = unconstrain_fixed unfix = unconstrain_fixed
def unconstrain_bounded(self, lower, upper): def unconstrain_bounded(self, lower, upper):
[param.unconstrain_bounded(lower, upper) for param in self.params] [param.unconstrain_bounded(lower, upper) for param in self.params]
unconstrain_bounded.__doc__ = Param.unconstrain_bounded.__doc__ unconstrain_bounded.__doc__ = Param.unconstrain_bounded.__doc__
def untie(self, *ties): def untie(self, *ties):
[param.untie(*ties) for param in self.params] [param.untie(*ties) for param in self.params]
__lt__ = lambda self, val: self._vals() < val __lt__ = lambda self, val: self._vals() < val
@ -547,11 +564,11 @@ class ParamConcatenation(object):
lx = max([p._max_len_values() for p in params]) lx = max([p._max_len_values() for p in params])
li = max([p._max_len_index(i) for p, i in itertools.izip(params, indices)]) li = max([p._max_len_index(i) for p, i in itertools.izip(params, indices)])
lt = max([p._max_len_names(tm, __tie_name__) for p, tm in itertools.izip(params, ties_matrices)]) lt = max([p._max_len_names(tm, __tie_name__) for p, tm in itertools.izip(params, ties_matrices)])
strings = [p.__str__(cm, i, tm, lc, lx, li, lt) for p, cm, i, tm in itertools.izip(params, constr_matrices, indices, ties_matrices)] strings = [p.__str__(cm, i, tm, lc, lx, li, lt) for p, cm, i, tm in itertools.izip(params,constr_matrices,indices,ties_matrices)]
return "\n".join(strings) return "\n".join(strings)
return "\n{}\n".format(" -" + "- | -".join(['-' * l for l in [li, lx, lc, lt]])).join(strings) return "\n{}\n".format(" -"+"- | -".join(['-'*l for l in [li,lx,lc,lt]])).join(strings)
def __repr__(self): def __repr__(self):
return "\n".join(map(repr, self.params)) return "\n".join(map(repr,self.params))
if __name__ == '__main__': if __name__ == '__main__':
@ -559,8 +576,8 @@ if __name__ == '__main__':
from GPy.core.parameterized import Parameterized from GPy.core.parameterized import Parameterized
from GPy.core.parameter import Param from GPy.core.parameter import Param
# X = numpy.random.randn(2,3,1,5,2,4,3) #X = numpy.random.randn(2,3,1,5,2,4,3)
X = numpy.random.randn(3, 2) X = numpy.random.randn(3,2)
print "random done" print "random done"
p = Param("q_mean", X) p = Param("q_mean", X)
p1 = Param("q_variance", numpy.random.rand(*p.shape)) p1 = Param("q_variance", numpy.random.rand(*p.shape))
@ -572,19 +589,19 @@ if __name__ == '__main__':
m = Parameterized() m = Parameterized()
rbf = Parameterized(name='rbf') rbf = Parameterized(name='rbf')
rbf.add_parameter(p3, p4) rbf.add_parameter(p3,p4)
m.add_parameter(p, p1, rbf) m.add_parameter(p,p1,rbf)
print "setting params" print "setting params"
# print m.q_v[3:5,[1,4,5]] #print m.q_v[3:5,[1,4,5]]
print "constraining variance" print "constraining variance"
# m[".*variance"].constrain_positive() #m[".*variance"].constrain_positive()
# print "constraining rbf" #print "constraining rbf"
# m.rbf_l.constrain_positive() #m.rbf_l.constrain_positive()
# m.q_variance[1,[0,5,11,19,2]].tie_to(m.rbf_v) #m.q_variance[1,[0,5,11,19,2]].tie_to(m.rbf_v)
# m.rbf_v.tie_to(m.rbf_l[0]) #m.rbf_v.tie_to(m.rbf_l[0])
# m.rbf_l[0].tie_to(m.rbf_l[1]) #m.rbf_l[0].tie_to(m.rbf_l[1])
# m.q_v.tie_to(m.rbf_v) #m.q_v.tie_to(m.rbf_v)
# m.rbf_l.tie_to(m.rbf_va) # m.rbf_l.tie_to(m.rbf_va)
# pt = numpy.array(params._get_params_transformed()) # pt = numpy.array(params._get_params_transformed())
# ptr = numpy.random.randn(*pt.shape) # ptr = numpy.random.randn(*pt.shape)

13
GPy/core/parameterization/parameter_core.py
View file

@ -52,6 +52,7 @@ class Parentable(object):
super(Parentable,self).__init__() super(Parentable,self).__init__()
self._direct_parent_ = direct_parent self._direct_parent_ = direct_parent
self._parent_index_ = parent_index self._parent_index_ = parent_index
self._highest_parent_ = highest_parent
def has_parent(self): def has_parent(self):
return self._direct_parent_ is not None return self._direct_parent_ is not None
@ -103,30 +104,30 @@ class Constrainable(Nameable):
if update: if update:
self.parameters_changed() self.parameters_changed()
def constrain_positive(self, warning=True): def constrain_positive(self, warning=True, update=True):
""" """
:param warning: print a warning if re-constraining parameters. :param warning: print a warning if re-constraining parameters.
Constrain this parameter to the default positive constraint. Constrain this parameter to the default positive constraint.
""" """
self.constrain(Logexp(), warning) self.constrain(Logexp(), warning=warning, update=update)
def constrain_negative(self, warning=True): def constrain_negative(self, warning=True, update=True):
""" """
:param warning: print a warning if re-constraining parameters. :param warning: print a warning if re-constraining parameters.
Constrain this parameter to the default negative constraint. Constrain this parameter to the default negative constraint.
""" """
self.constrain(NegativeLogexp(), warning) self.constrain(NegativeLogexp(), warning=warning, update=update)
def constrain_bounded(self, lower, upper, warning=True): def constrain_bounded(self, lower, upper, warning=True, update=True):
""" """
:param lower, upper: the limits to bound this parameter to :param lower, upper: the limits to bound this parameter to
:param warning: print a warning if re-constraining parameters. :param warning: print a warning if re-constraining parameters.
Constrain this parameter to lie within the given range. Constrain this parameter to lie within the given range.
""" """
self.constrain(Logistic(lower, upper), warning) self.constrain(Logistic(lower, upper), warning=warning, update=update)
def unconstrain(self, *transforms): def unconstrain(self, *transforms):
""" """

15
GPy/core/parameterization/parameterized.py
View file

@ -248,6 +248,16 @@ class Parameterized(Constrainable, Pickleable, Observable):
cPickle.dump(self, f, protocol) cPickle.dump(self, f, protocol)
def copy(self): def copy(self):
"""Returns a (deep) copy of the current model """ """Returns a (deep) copy of the current model """
#dc = dict()
#for k, v in self.__dict__.iteritems():
#if k not in ['_highest_parent_', '_direct_parent_']:
#dc[k] = copy.deepcopy(v)
#dc = copy.deepcopy(self.__dict__)
#dc['_highest_parent_'] = None
#dc['_direct_parent_'] = None
#s = self.__class__.new()
#s.__dict__ = dc
return copy.deepcopy(self) return copy.deepcopy(self)
def __getstate__(self): def __getstate__(self):
if self._has_get_set_state(): if self._has_get_set_state():
@ -413,6 +423,8 @@ class Parameterized(Constrainable, Pickleable, Observable):
#=========================================================================== #===========================================================================
# Convenience for fixed, tied checking of param: # Convenience for fixed, tied checking of param:
#=========================================================================== #===========================================================================
def fixed_indices(self):
return np.array([x.is_fixed for x in self._parameters_])
def _is_fixed(self, param): def _is_fixed(self, param):
# returns if the whole param is fixed # returns if the whole param is fixed
if not self._has_fixes(): if not self._has_fixes():
@ -442,7 +454,8 @@ class Parameterized(Constrainable, Pickleable, Observable):
# if removing constraints before adding new is not wanted, just delete the above line! # if removing constraints before adding new is not wanted, just delete the above line!
self.constraints.add(transform, rav_i) self.constraints.add(transform, rav_i)
param = self._get_original(param) param = self._get_original(param)
param._set_params(transform.initialize(param._get_params())) if not (transform == __fixed__):
param._set_params(transform.initialize(param._get_params()), update=False)
if warning and any(reconstrained): if warning and any(reconstrained):
# if you want to print the whole params object, which was reconstrained use: # if you want to print the whole params object, which was reconstrained use:
# m = str(param[self._backtranslate_index(param, reconstrained)]) # m = str(param[self._backtranslate_index(param, reconstrained)])

4
GPy/core/parameterization/transformations.py
View file

@ -162,7 +162,9 @@ class Logistic(Transformation):
def initialize(self, f): def initialize(self, f):
if np.any(np.logical_or(f < self.lower, f > self.upper)): if np.any(np.logical_or(f < self.lower, f > self.upper)):
print "Warning: changing parameters to satisfy constraints" print "Warning: changing parameters to satisfy constraints"
return np.where(np.logical_or(f < self.lower, f > self.upper), self.f(f * 0.), f) #return np.where(np.logical_or(f < self.lower, f > self.upper), self.f(f * 0.), f)
#FIXME: Max, zeros_like right?
return np.where(np.logical_or(f < self.lower, f > self.upper), self.f(np.zeros_like(f)), f)
def __str__(self): def __str__(self):
return '{},{}'.format(self.lower, self.upper) return '{},{}'.format(self.lower, self.upper)

100
GPy/inference/latent_function_inference/laplace.py
View file

@ -32,7 +32,7 @@ class LaplaceInference(object):
self._mode_finding_tolerance = 1e-7 self._mode_finding_tolerance = 1e-7
self._mode_finding_max_iter = 40 self._mode_finding_max_iter = 40
self.bad_fhat = True self.bad_fhat = True
self._previous_Ki_fhat = None
def inference(self, kern, X, likelihood, Y, Y_metadata=None): def inference(self, kern, X, likelihood, Y, Y_metadata=None):
""" """
@ -50,16 +50,17 @@ class LaplaceInference(object):
Ki_f_init = np.zeros_like(Y) Ki_f_init = np.zeros_like(Y)
else: else:
Ki_f_init = self._previous_Ki_fhat Ki_f_init = self._previous_Ki_fhat
f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata) f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)
#Compute hessian and other variables at mode #Compute hessian and other variables at mode
log_marginal, Ki_W_i, K_Wi_i, dL_dK, woodbury_vector = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, Y_metadata) log_marginal, woodbury_vector, woodbury_inv, dL_dK, dL_dthetaL = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, kern, Y_metadata)
#likelihood.gradient = self.likelihood_gradients()
kern.update_gradients_full(dL_dK, X) kern.update_gradients_full(dL_dK, X)
likelihood.update_gradients(dL_dthetaL)
self._previous_Ki_fhat = Ki_fhat.copy() self._previous_Ki_fhat = Ki_fhat.copy()
return Posterior(woodbury_vector=woodbury_vector, woodbury_inv = K_Wi_i, K=K), log_marginal, {'dL_dK':dL_dK} return Posterior(woodbury_vector=woodbury_vector, woodbury_inv=woodbury_inv, K=K), log_marginal, {'dL_dK':dL_dK}
def rasm_mode(self, K, Y, likelihood, Ki_f_init, Y_metadata=None): def rasm_mode(self, K, Y, likelihood, Ki_f_init, Y_metadata=None):
""" """
@ -133,13 +134,15 @@ class LaplaceInference(object):
return f, Ki_f return f, Ki_f
def mode_computations(self, f_hat, Ki_f, K, Y, likelihood, kern, Y_metadata):
def mode_computations(self, f_hat, Ki_f, K, Y, likelihood, Y_metadata):
""" """
At the mode, compute the hessian and effective covariance matrix. At the mode, compute the hessian and effective covariance matrix.
returns: logZ : approximation to the marginal likelihood returns: logZ : approximation to the marginal likelihood
Cov : the approximation to the covariance matrix woodbury_vector : variable required for calculating the approximation to the covariance matrix
woodbury_inv : variable required for calculating the approximation to the covariance matrix
dL_dthetaL : array of derivatives (1 x num_kernel_params)
dL_dthetaL : array of derivatives (1 x num_likelihood_params)
""" """
#At this point get the hessian matrix (or vector as W is diagonal) #At this point get the hessian matrix (or vector as W is diagonal)
W = -likelihood.d2logpdf_df2(f_hat, Y, extra_data=Y_metadata) W = -likelihood.d2logpdf_df2(f_hat, Y, extra_data=Y_metadata)
@ -153,45 +156,62 @@ class LaplaceInference(object):
#compute the log marginal #compute the log marginal
log_marginal = -0.5*np.dot(Ki_f.flatten(), f_hat.flatten()) + likelihood.logpdf(f_hat, Y, extra_data=Y_metadata) - np.sum(np.log(np.diag(L))) log_marginal = -0.5*np.dot(Ki_f.flatten(), f_hat.flatten()) + likelihood.logpdf(f_hat, Y, extra_data=Y_metadata) - np.sum(np.log(np.diag(L)))
#compute dL_dK #Compute vival matrices for derivatives
explicit_part = 0.5*(np.dot(Ki_f, Ki_f.T) - K_Wi_i) dW_df = -likelihood.d3logpdf_df3(f_hat, Y, extra_data=Y_metadata) # -d3lik_d3fhat
#Implicit
d3lik_d3fhat = likelihood.d3logpdf_df3(f_hat, Y, extra_data=Y_metadata)
dL_dfhat = 0.5*(np.diag(Ki_W_i)[:, None]*d3lik_d3fhat) #why isn't this -0.5? s2 in R&W p126 line 9.
woodbury_vector = likelihood.dlogpdf_df(f_hat, Y, extra_data=Y_metadata) woodbury_vector = likelihood.dlogpdf_df(f_hat, Y, extra_data=Y_metadata)
implicit_part = np.dot(woodbury_vector, dL_dfhat.T).dot(np.eye(Y.shape[0]) - np.dot(K, K_Wi_i)) dL_dfhat = -0.5*(np.diag(Ki_W_i)[:, None]*dW_df) #why isn't this -0.5? s2 in R&W p126 line 9.
#BiK, _ = dpotrs(L, K, lower=1)
#dL_dfhat = 0.5*np.diag(BiK)[:, None]*dW_df
I_KW_i = np.eye(Y.shape[0]) - np.dot(K, K_Wi_i)
dL_dK = explicit_part + implicit_part ####################
#compute dL_dK#
return log_marginal, Ki_W_i, K_Wi_i, dL_dK, woodbury_vector ####################
if kern.size > 0 and not kern.is_fixed:
def likelihood_gradients(self):
"""
Gradients with respect to likelihood parameters (dL_dthetaL)
:rtype: array of derivatives (1 x num_likelihood_params)
"""
dL_dfhat, I_KW_i = self._shared_gradients_components()
dlik_dthetaL, dlik_grad_dthetaL, dlik_hess_dthetaL = likelihood._laplace_gradients(self.f_hat, self.data, extra_data=self.extra_data)
num_params = len(self._get_param_names())
# make space for one derivative for each likelihood parameter
dL_dthetaL = np.zeros(num_params)
for thetaL_i in range(num_params):
#Explicit #Explicit
dL_dthetaL_exp = ( np.sum(dlik_dthetaL[:, thetaL_i]) explicit_part = 0.5*(np.dot(Ki_f, Ki_f.T) - K_Wi_i)
#- 0.5*np.trace(mdot(self.Ki_W_i, (self.K, np.diagflat(dlik_hess_dthetaL[thetaL_i]))))
+ np.dot(0.5*np.diag(self.Ki_W_i)[:,None].T, dlik_hess_dthetaL[:, thetaL_i])
)
#Implicit #Implicit
dfhat_dthetaL = mdot(I_KW_i, self.K, dlik_grad_dthetaL[:, thetaL_i]) implicit_part = np.dot(woodbury_vector, dL_dfhat.T).dot(I_KW_i)
dL_dthetaL_imp = np.dot(dL_dfhat, dfhat_dthetaL)
dL_dthetaL[thetaL_i] = dL_dthetaL_exp + dL_dthetaL_imp
return dL_dthetaL dL_dK = explicit_part + implicit_part
else:
dL_dK = np.zeros(likelihood.size)
####################
#compute dL_dthetaL#
####################
if likelihood.size > 0 and not likelihood.is_fixed:
dlik_dthetaL, dlik_grad_dthetaL, dlik_hess_dthetaL = likelihood._laplace_gradients(f_hat, Y, extra_data=Y_metadata)
num_params = likelihood.size
# make space for one derivative for each likelihood parameter
dL_dthetaL = np.zeros(num_params)
for thetaL_i in range(num_params):
#Explicit
dL_dthetaL_exp = ( np.sum(dlik_dthetaL[thetaL_i])
# The + comes from the fact that dlik_hess_dthetaL == -dW_dthetaL
+ 0.5*np.sum(np.diag(Ki_W_i).flatten()*dlik_hess_dthetaL[:, thetaL_i].flatten())
)
#Implicit
dfhat_dthetaL = mdot(I_KW_i, K, dlik_grad_dthetaL[:, thetaL_i])
#dfhat_dthetaL = mdot(Ki_W_i, dlik_grad_dthetaL[:, thetaL_i])
dL_dthetaL_imp = np.dot(dL_dfhat.T, dfhat_dthetaL)
dL_dthetaL[thetaL_i] = dL_dthetaL_exp + dL_dthetaL_imp
else:
dL_dthetaL = np.zeros(likelihood.size)
return log_marginal, woodbury_vector, K_Wi_i, dL_dK, dL_dthetaL
#def likelihood_gradients(self, f_hat, K, Y, Ki_W_i, dL_dfhat, I_KW_i, likelihood, Y_metadata):
#"""
#Gradients with respect to likelihood parameters (dL_dthetaL)
#:rtype: array of derivatives (1 x num_likelihood_params)
#"""
def _compute_B_statistics(self, K, W, log_concave): def _compute_B_statistics(self, K, W, log_concave):
""" """

17
GPy/likelihoods/gaussian.py
View file

@ -130,7 +130,10 @@ class Gaussian(Likelihood):
:rtype: float :rtype: float
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
return -0.5*(np.sum((y-link_f)**2/self.variance) + self.ln_det_K + self.N*np.log(2.*np.pi)) N = y.shape[0]
ln_det_cov = N*np.log(self.variance)
return -0.5*(np.sum((y-link_f)**2/self.variance) + ln_det_cov + N*np.log(2.*np.pi))
def dlogpdf_dlink(self, link_f, y, extra_data=None): def dlogpdf_dlink(self, link_f, y, extra_data=None):
""" """
@ -175,7 +178,8 @@ class Gaussian(Likelihood):
(the distribution for y_i depends only on link(f_i) not on link(f_(j!=i)) (the distribution for y_i depends only on link(f_i) not on link(f_(j!=i))
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
hess = -(1.0/self.variance)*np.ones((self.N, 1)) N = y.shape[0]
hess = -(1.0/self.variance)*np.ones((N, 1))
return hess return hess
def d3logpdf_dlink3(self, link_f, y, extra_data=None): def d3logpdf_dlink3(self, link_f, y, extra_data=None):
@ -194,7 +198,8 @@ class Gaussian(Likelihood):
:rtype: Nx1 array :rtype: Nx1 array
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
d3logpdf_dlink3 = np.diagonal(0*self.I)[:, None] N = y.shape[0]
d3logpdf_dlink3 = np.zeros((N,1))
return d3logpdf_dlink3 return d3logpdf_dlink3
def dlogpdf_link_dvar(self, link_f, y, extra_data=None): def dlogpdf_link_dvar(self, link_f, y, extra_data=None):
@ -215,7 +220,8 @@ class Gaussian(Likelihood):
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
e = y - link_f e = y - link_f
s_4 = 1.0/(self.variance**2) s_4 = 1.0/(self.variance**2)
dlik_dsigma = -0.5*self.N/self.variance + 0.5*s_4*np.sum(np.square(e)) N = y.shape[0]
dlik_dsigma = -0.5*N/self.variance + 0.5*s_4*np.sum(np.square(e))
return np.sum(dlik_dsigma) # Sure about this sum? return np.sum(dlik_dsigma) # Sure about this sum?
def dlogpdf_dlink_dvar(self, link_f, y, extra_data=None): def dlogpdf_dlink_dvar(self, link_f, y, extra_data=None):
@ -255,7 +261,8 @@ class Gaussian(Likelihood):
""" """
assert np.asarray(link_f).shape == np.asarray(y).shape assert np.asarray(link_f).shape == np.asarray(y).shape
s_4 = 1.0/(self.variance**2) s_4 = 1.0/(self.variance**2)
d2logpdf_dlink2_dvar = np.diag(s_4*self.I)[:, None] N = y.shape[0]
d2logpdf_dlink2_dvar = np.ones((N,1))*s_4
return d2logpdf_dlink2_dvar return d2logpdf_dlink2_dvar
def dlogpdf_link_dtheta(self, f, y, extra_data=None): def dlogpdf_link_dtheta(self, f, y, extra_data=None):

25
GPy/likelihoods/likelihood.py
View file

@ -44,6 +44,10 @@ class Likelihood(Parameterized):
def _gradients(self,partial): def _gradients(self,partial):
return np.zeros(0) return np.zeros(0)
def update_gradients(self, partial):
if self.size > 0:
raise NotImplementedError('Must be implemented for likelihoods with parameters to be optimized')
def _preprocess_values(self,Y): def _preprocess_values(self,Y):
""" """
In case it is needed, this function assess the output values or makes any pertinent transformation on them. In case it is needed, this function assess the output values or makes any pertinent transformation on them.
@ -303,31 +307,31 @@ class Likelihood(Parameterized):
""" """
TODO: Doc strings TODO: Doc strings
""" """
if len(self._get_param_names()) > 0: if self.size > 0:
link_f = self.gp_link.transf(f) link_f = self.gp_link.transf(f)
return self.dlogpdf_link_dtheta(link_f, y, extra_data=extra_data) return self.dlogpdf_link_dtheta(link_f, y, extra_data=extra_data)
else: else:
#Is no parameters so return an empty array for its derivatives #Is no parameters so return an empty array for its derivatives
return np.empty([1, 0]) return np.zeros([1, 0])
def dlogpdf_df_dtheta(self, f, y, extra_data=None): def dlogpdf_df_dtheta(self, f, y, extra_data=None):
""" """
TODO: Doc strings TODO: Doc strings
""" """
if len(self._get_param_names()) > 0: if self.size > 0:
link_f = self.gp_link.transf(f) link_f = self.gp_link.transf(f)
dlink_df = self.gp_link.dtransf_df(f) dlink_df = self.gp_link.dtransf_df(f)
dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(link_f, y, extra_data=extra_data) dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(link_f, y, extra_data=extra_data)
return chain_1(dlogpdf_dlink_dtheta, dlink_df) return chain_1(dlogpdf_dlink_dtheta, dlink_df)
else: else:
#Is no parameters so return an empty array for its derivatives #Is no parameters so return an empty array for its derivatives
return np.empty([f.shape[0], 0]) return np.zeros([f.shape[0], 0])
def d2logpdf_df2_dtheta(self, f, y, extra_data=None): def d2logpdf_df2_dtheta(self, f, y, extra_data=None):
""" """
TODO: Doc strings TODO: Doc strings
""" """
if len(self._get_param_names()) > 0: if self.size > 0:
link_f = self.gp_link.transf(f) link_f = self.gp_link.transf(f)
dlink_df = self.gp_link.dtransf_df(f) dlink_df = self.gp_link.dtransf_df(f)
d2link_df2 = self.gp_link.d2transf_df2(f) d2link_df2 = self.gp_link.d2transf_df2(f)
@ -336,7 +340,7 @@ class Likelihood(Parameterized):
return chain_2(d2logpdf_dlink2_dtheta, dlink_df, dlogpdf_dlink_dtheta, d2link_df2) return chain_2(d2logpdf_dlink2_dtheta, dlink_df, dlogpdf_dlink_dtheta, d2link_df2)
else: else:
#Is no parameters so return an empty array for its derivatives #Is no parameters so return an empty array for its derivatives
return np.empty([f.shape[0], 0]) return np.zeros([f.shape[0], 0])
def _laplace_gradients(self, f, y, extra_data=None): def _laplace_gradients(self, f, y, extra_data=None):
dlogpdf_dtheta = self.dlogpdf_dtheta(f, y, extra_data=extra_data) dlogpdf_dtheta = self.dlogpdf_dtheta(f, y, extra_data=extra_data)
@ -345,9 +349,12 @@ class Likelihood(Parameterized):
#Parameters are stacked vertically. Must be listed in same order as 'get_param_names' #Parameters are stacked vertically. Must be listed in same order as 'get_param_names'
# ensure we have gradients for every parameter we want to optimize # ensure we have gradients for every parameter we want to optimize
assert dlogpdf_dtheta.shape[1] == len(self._get_param_names()) try:
assert dlogpdf_df_dtheta.shape[1] == len(self._get_param_names()) assert len(dlogpdf_dtheta) == self.size #1 x num_param array
assert d2logpdf_df2_dtheta.shape[1] == len(self._get_param_names()) assert dlogpdf_df_dtheta.shape[1] == self.size #f x num_param matrix
assert d2logpdf_df2_dtheta.shape[1] == self.size #f x num_param matrix
except Exception as e:
import ipdb; ipdb.set_trace() # XXX BREAKPOINT
return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta

56
GPy/likelihoods/student_t.py
View file

@ -8,6 +8,7 @@ import link_functions
from scipy import stats, integrate from scipy import stats, integrate
from scipy.special import gammaln, gamma from scipy.special import gammaln, gamma
from likelihood import Likelihood from likelihood import Likelihood
from ..core.parameterization import Param
class StudentT(Likelihood): class StudentT(Likelihood):
""" """
@ -19,26 +20,30 @@ class StudentT(Likelihood):
p(y_{i}|\\lambda(f_{i})) = \\frac{\\Gamma\\left(\\frac{v+1}{2}\\right)}{\\Gamma\\left(\\frac{v}{2}\\right)\\sqrt{v\\pi\\sigma^{2}}}\\left(1 + \\frac{1}{v}\\left(\\frac{(y_{i} - f_{i})^{2}}{\\sigma^{2}}\\right)\\right)^{\\frac{-v+1}{2}} p(y_{i}|\\lambda(f_{i})) = \\frac{\\Gamma\\left(\\frac{v+1}{2}\\right)}{\\Gamma\\left(\\frac{v}{2}\\right)\\sqrt{v\\pi\\sigma^{2}}}\\left(1 + \\frac{1}{v}\\left(\\frac{(y_{i} - f_{i})^{2}}{\\sigma^{2}}\\right)\\right)^{\\frac{-v+1}{2}}
""" """
def __init__(self,gp_link=None,analytical_mean=True,analytical_variance=True, deg_free=5, sigma2=2): def __init__(self,gp_link=None, deg_free=5, sigma2=2):
self.v = deg_free if gp_link is None:
self.sigma2 = sigma2 gp_link = link_functions.Identity()
super(StudentT, self).__init__(gp_link, name='Student_T')
self.sigma2 = Param('t_noise', float(sigma2))
self.v = Param('deg_free', float(deg_free))
self.add_parameter(self.sigma2)
self.add_parameter(self.v)
self.v.constrain_fixed()
self._set_params(np.asarray(sigma2))
super(StudentT, self).__init__(gp_link,analytical_mean,analytical_variance)
self.log_concave = False self.log_concave = False
def _get_params(self): def parameters_changed(self):
return np.asarray(self.sigma2) self.variance = (self.v / float(self.v - 2)) * self.sigma2
def _get_param_names(self): def update_gradients(self, partial):
return ["t_noise_std2"] """
Pull out the gradients, be careful as the order must match the order
def _set_params(self, x): in which the parameters are added
self.sigma2 = float(x) """
self.sigma2.gradient = partial[0]
@property self.v.gradient = partial[1]
def variance(self, extra_data=None):
return (self.v / float(self.v - 2)) * self.sigma2
def pdf_link(self, link_f, y, extra_data=None): def pdf_link(self, link_f, y, extra_data=None):
""" """
@ -82,10 +87,14 @@ class StudentT(Likelihood):
""" """
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
e = y - link_f e = y - link_f
#FIXME:
#Why does np.log(1 + (1/self.v)*((y-link_f)**2)/self.sigma2) suppress the divide by zero?!
#But np.log(1 + (1/float(self.v))*((y-link_f)**2)/self.sigma2) throws it correctly
#print - 0.5*(self.v + 1)*np.log(1 + (1/np.float(self.v))*((e**2)/self.sigma2))
objective = (+ gammaln((self.v + 1) * 0.5) objective = (+ gammaln((self.v + 1) * 0.5)
- gammaln(self.v * 0.5) - gammaln(self.v * 0.5)
- 0.5*np.log(self.sigma2 * self.v * np.pi) - 0.5*np.log(self.sigma2 * self.v * np.pi)
- 0.5*(self.v + 1)*np.log(1 + (1/np.float(self.v))*((e**2)/self.sigma2)) - 0.5*(self.v + 1)*np.log(1 + (1/np.float(self.v))*((e**2)/self.sigma2))
) )
return np.sum(objective) return np.sum(objective)
@ -222,15 +231,18 @@ class StudentT(Likelihood):
def dlogpdf_link_dtheta(self, f, y, extra_data=None): def dlogpdf_link_dtheta(self, f, y, extra_data=None):
dlogpdf_dvar = self.dlogpdf_link_dvar(f, y, extra_data=extra_data) dlogpdf_dvar = self.dlogpdf_link_dvar(f, y, extra_data=extra_data)
return np.asarray([[dlogpdf_dvar]]) dlogpdf_dv = np.zeros_like(dlogpdf_dvar) #FIXME: Not done yet
return np.hstack((dlogpdf_dvar, dlogpdf_dv))
def dlogpdf_dlink_dtheta(self, f, y, extra_data=None): def dlogpdf_dlink_dtheta(self, f, y, extra_data=None):
dlogpdf_dlink_dvar = self.dlogpdf_dlink_dvar(f, y, extra_data=extra_data) dlogpdf_dlink_dvar = self.dlogpdf_dlink_dvar(f, y, extra_data=extra_data)
return dlogpdf_dlink_dvar dlogpdf_dlink_dv = np.zeros_like(dlogpdf_dlink_dvar) #FIXME: Not done yet
return np.hstack((dlogpdf_dlink_dvar, dlogpdf_dlink_dv))
def d2logpdf_dlink2_dtheta(self, f, y, extra_data=None): def d2logpdf_dlink2_dtheta(self, f, y, extra_data=None):
d2logpdf_dlink2_dvar = self.d2logpdf_dlink2_dvar(f, y, extra_data=extra_data) d2logpdf_dlink2_dvar = self.d2logpdf_dlink2_dvar(f, y, extra_data=extra_data)
return d2logpdf_dlink2_dvar d2logpdf_dlink2_dv = np.zeros_like(d2logpdf_dlink2_dvar) #FIXME: Not done yet
return np.hstack((d2logpdf_dlink2_dvar, d2logpdf_dlink2_dv))
def _predictive_variance_analytical(self, mu, sigma, predictive_mean=None): def _predictive_variance_analytical(self, mu, sigma, predictive_mean=None):
""" """

37
GPy/models/gradient_checker.py
View file

@ -1,9 +1,10 @@
# ## Copyright (c) 2012, GPy authors (see AUTHORS.txt). # ## Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt) # Licensed under the BSD 3-clause license (see LICENSE.txt)
from GPy.core.model import Model from ..core.model import Model
import itertools import itertools
import numpy import numpy
from ..core.parameterization import Param
def get_shape(x): def get_shape(x):
if isinstance(x, numpy.ndarray): if isinstance(x, numpy.ndarray):
@ -59,7 +60,7 @@ class GradientChecker(Model):
grad.randomize() grad.randomize()
grad.checkgrad(verbose=1) grad.checkgrad(verbose=1)
""" """
Model.__init__(self) Model.__init__(self, 'GradientChecker')
if isinstance(x0, (list, tuple)) and names is None: if isinstance(x0, (list, tuple)) and names is None:
self.shapes = [get_shape(xi) for xi in x0] self.shapes = [get_shape(xi) for xi in x0]
self.names = ['X{i}'.format(i=i) for i in range(len(x0))] self.names = ['X{i}'.format(i=i) for i in range(len(x0))]
@ -72,8 +73,10 @@ class GradientChecker(Model):
else: else:
self.names = names self.names = names
self.shapes = [get_shape(x0)] self.shapes = [get_shape(x0)]
for name, xi in zip(self.names, at_least_one_element(x0)): for name, xi in zip(self.names, at_least_one_element(x0)):
self.__setattr__(name, xi) self.__setattr__(name, Param(name, xi))
self.add_parameter(self.__getattribute__(name))
# self._param_names = [] # self._param_names = []
# for name, shape in zip(self.names, self.shapes): # for name, shape in zip(self.names, self.shapes):
# self._param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape)))))) # self._param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape))))))
@ -93,20 +96,18 @@ class GradientChecker(Model):
def _log_likelihood_gradients(self): def _log_likelihood_gradients(self):
return numpy.atleast_1d(self.df(*self._get_x(), **self.kwargs)).flatten() return numpy.atleast_1d(self.df(*self._get_x(), **self.kwargs)).flatten()
#def _get_params(self):
#return numpy.atleast_1d(numpy.hstack(map(lambda name: flatten_if_needed(self.__getattribute__(name)), self.names)))
def _get_params(self): #def _set_params(self, x):
return numpy.atleast_1d(numpy.hstack(map(lambda name: flatten_if_needed(self.__getattribute__(name)), self.names))) #current_index = 0
#for name, shape in zip(self.names, self.shapes):
#current_size = numpy.prod(shape)
#self.__setattr__(name, x[current_index:current_index + current_size].reshape(shape))
#current_index += current_size
#def _get_param_names(self):
def _set_params(self, x): #_param_names = []
current_index = 0 #for name, shape in zip(self.names, self.shapes):
for name, shape in zip(self.names, self.shapes): #_param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape))))))
current_size = numpy.prod(shape) #return _param_names
self.__setattr__(name, x[current_index:current_index + current_size].reshape(shape))
current_index += current_size
def _get_param_names(self):
_param_names = []
for name, shape in zip(self.names, self.shapes):
_param_names.extend(map(lambda nameshape: ('_'.join(nameshape)).strip('_'), itertools.izip(itertools.repeat(name), itertools.imap(lambda t: '_'.join(map(str, t)), itertools.product(*map(lambda xi: range(xi), shape))))))
return _param_names

235
GPy/testing/likelihood_tests.py
View file

@ -4,10 +4,11 @@ import GPy
from GPy.models import GradientChecker from GPy.models import GradientChecker
import functools import functools
import inspect import inspect
from GPy.likelihoods.noise_models import gp_transformations from GPy.likelihoods import link_functions
from ..core.parameterization import Param
from functools import partial from functools import partial
#np.random.seed(300) #np.random.seed(300)
np.random.seed(7) #np.random.seed(7)
def dparam_partial(inst_func, *args): def dparam_partial(inst_func, *args):
""" """
@ -22,12 +23,14 @@ def dparam_partial(inst_func, *args):
the f or Y that are being used in the function whilst we tweak the the f or Y that are being used in the function whilst we tweak the
param param
""" """
def param_func(param, inst_func, args): def param_func(param_val, param_name, inst_func, args):
inst_func.im_self._set_params(param) #inst_func.im_self._set_params(param)
#inst_func.im_self.add_parameter(Param(param_name, param_val))
inst_func.im_self[param_name] = param_val
return inst_func(*args) return inst_func(*args)
return functools.partial(param_func, inst_func=inst_func, args=args) return functools.partial(param_func, inst_func=inst_func, args=args)
def dparam_checkgrad(func, dfunc, params, args, constraints=None, randomize=False, verbose=False): def dparam_checkgrad(func, dfunc, params, params_names, args, constraints=None, randomize=False, verbose=False):
""" """
checkgrad expects a f: R^N -> R^1 and df: R^N -> R^N checkgrad expects a f: R^N -> R^1 and df: R^N -> R^N
However if we are holding other parameters fixed and moving something else However if we are holding other parameters fixed and moving something else
@ -38,27 +41,34 @@ def dparam_checkgrad(func, dfunc, params, args, constraints=None, randomize=Fals
The number of parameters and N is the number of data The number of parameters and N is the number of data
Need to take a slice out from f and a slice out of df Need to take a slice out from f and a slice out of df
""" """
#print "\n{} likelihood: {} vs {}".format(func.im_self.__class__.__name__, print "\n{} likelihood: {} vs {}".format(func.im_self.__class__.__name__,
#func.__name__, dfunc.__name__) func.__name__, dfunc.__name__)
partial_f = dparam_partial(func, *args) partial_f = dparam_partial(func, *args)
partial_df = dparam_partial(dfunc, *args) partial_df = dparam_partial(dfunc, *args)
gradchecking = True gradchecking = True
for param in params: zipped_params = zip(params, params_names)
fnum = np.atleast_1d(partial_f(param)).shape[0] for param_ind, (param_val, param_name) in enumerate(zipped_params):
dfnum = np.atleast_1d(partial_df(param)).shape[0] #Check one parameter at a time, make sure it is 2d (as some gradients only return arrays) then strip out the parameter
fnum = np.atleast_2d(partial_f(param_val, param_name))[:, param_ind].shape[0]
dfnum = np.atleast_2d(partial_df(param_val, param_name))[:, param_ind].shape[0]
for fixed_val in range(dfnum): for fixed_val in range(dfnum):
#dlik and dlik_dvar gives back 1 value for each #dlik and dlik_dvar gives back 1 value for each
f_ind = min(fnum, fixed_val+1) - 1 f_ind = min(fnum, fixed_val+1) - 1
print "fnum: {} dfnum: {} f_ind: {} fixed_val: {}".format(fnum, dfnum, f_ind, fixed_val) print "fnum: {} dfnum: {} f_ind: {} fixed_val: {}".format(fnum, dfnum, f_ind, fixed_val)
#Make grad checker with this param moving, note that set_params is NOT being called #Make grad checker with this param moving, note that set_params is NOT being called
#The parameter is being set directly with __setattr__ #The parameter is being set directly with __setattr__
grad = GradientChecker(lambda x: np.atleast_1d(partial_f(x))[f_ind], #Check only the parameter and function value we wish to check at a time
lambda x : np.atleast_1d(partial_df(x))[fixed_val], grad = GradientChecker(lambda p_val: np.atleast_2d(partial_f(p_val, param_name))[f_ind, param_ind],
param, 'p') lambda p_val: np.atleast_2d(partial_df(p_val, param_name))[fixed_val, param_ind],
#This is not general for more than one param... param_val, [param_name])
if constraints is not None: if constraints is not None:
for constraint in constraints: for constrain_param, constraint in constraints:
constraint('p', grad) if grad.grep_param_names(constrain_param):
constraint(constrain_param, grad)
else:
print "parameter didn't exist"
print constrain_param, " ", constraint
if randomize: if randomize:
grad.randomize() grad.randomize()
if verbose: if verbose:
@ -107,17 +117,20 @@ class TestNoiseModels(object):
#################################################### ####################################################
# Constraint wrappers so we can just list them off # # Constraint wrappers so we can just list them off #
#################################################### ####################################################
def constrain_fixed(regex, model):
model[regex].constrain_fixed()
def constrain_negative(regex, model): def constrain_negative(regex, model):
model.constrain_negative(regex) model[regex].constrain_negative()
def constrain_positive(regex, model): def constrain_positive(regex, model):
model.constrain_positive(regex) model[regex].constrain_positive()
def constrain_bounded(regex, model, lower, upper): def constrain_bounded(regex, model, lower, upper):
""" """
Used like: partial(constrain_bounded, lower=0, upper=1) Used like: partial(constrain_bounded, lower=0, upper=1)
""" """
model.constrain_bounded(regex, lower, upper) model[regex].constrain_bounded(lower, upper)
""" """
Dictionary where we nest models we would like to check Dictionary where we nest models we would like to check
@ -134,71 +147,72 @@ class TestNoiseModels(object):
} }
""" """
noise_models = {"Student_t_default": { noise_models = {"Student_t_default": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
#"constraints": [("t_noise", constrain_positive), ("deg_free", partial(constrain_fixed, value=5))]
}, },
"laplace": True "laplace": True
}, },
"Student_t_1_var": { "Student_t_1_var": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [1.0], "vals": [1.0],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_small_var": { "Student_t_small_var": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [0.01], "vals": [0.01],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_large_var": { "Student_t_large_var": {
"model": GPy.likelihoods.student_t(deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [10.0], "vals": [10.0],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_approx_gauss": { "Student_t_approx_gauss": {
"model": GPy.likelihoods.student_t(deg_free=1000, sigma2=self.var), "model": GPy.likelihoods.StudentT(deg_free=1000, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
}, },
"laplace": True "laplace": True
}, },
"Student_t_log": { "Student_t_log": {
"model": GPy.likelihoods.student_t(gp_link=gp_transformations.Log(), deg_free=5, sigma2=self.var), "model": GPy.likelihoods.StudentT(gp_link=link_functions.Log(), deg_free=5, sigma2=self.var),
"grad_params": { "grad_params": {
"names": ["t_noise"], "names": ["t_noise"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
}, },
"laplace": True "laplace": True
}, },
"Gaussian_default": { "Gaussian_default": {
"model": GPy.likelihoods.gaussian(variance=self.var, D=self.D, N=self.N), "model": GPy.likelihoods.Gaussian(variance=self.var),
"grad_params": { "grad_params": {
"names": ["noise_model_variance"], "names": ["variance"],
"vals": [self.var], "vals": [self.var],
"constraints": [constrain_positive] "constraints": [("variance", constrain_positive)]
}, },
"laplace": True, "laplace": True,
"ep": True "ep": True
}, },
#"Gaussian_log": { #"Gaussian_log": {
#"model": GPy.likelihoods.gaussian(gp_link=gp_transformations.Log(), variance=self.var, D=self.D, N=self.N), #"model": GPy.likelihoods.gaussian(gp_link=link_functions.Log(), variance=self.var, D=self.D, N=self.N),
#"grad_params": { #"grad_params": {
#"names": ["noise_model_variance"], #"names": ["noise_model_variance"],
#"vals": [self.var], #"vals": [self.var],
@ -207,7 +221,7 @@ class TestNoiseModels(object):
#"laplace": True #"laplace": True
#}, #},
#"Gaussian_probit": { #"Gaussian_probit": {
#"model": GPy.likelihoods.gaussian(gp_link=gp_transformations.Probit(), variance=self.var, D=self.D, N=self.N), #"model": GPy.likelihoods.gaussian(gp_link=link_functions.Probit(), variance=self.var, D=self.D, N=self.N),
#"grad_params": { #"grad_params": {
#"names": ["noise_model_variance"], #"names": ["noise_model_variance"],
#"vals": [self.var], #"vals": [self.var],
@ -216,7 +230,7 @@ class TestNoiseModels(object):
#"laplace": True #"laplace": True
#}, #},
#"Gaussian_log_ex": { #"Gaussian_log_ex": {
#"model": GPy.likelihoods.gaussian(gp_link=gp_transformations.Log_ex_1(), variance=self.var, D=self.D, N=self.N), #"model": GPy.likelihoods.gaussian(gp_link=link_functions.Log_ex_1(), variance=self.var, D=self.D, N=self.N),
#"grad_params": { #"grad_params": {
#"names": ["noise_model_variance"], #"names": ["noise_model_variance"],
#"vals": [self.var], #"vals": [self.var],
@ -225,31 +239,31 @@ class TestNoiseModels(object):
#"laplace": True #"laplace": True
#}, #},
"Bernoulli_default": { "Bernoulli_default": {
"model": GPy.likelihoods.bernoulli(), "model": GPy.likelihoods.Bernoulli(),
"link_f_constraints": [partial(constrain_bounded, lower=0, upper=1)], "link_f_constraints": [partial(constrain_bounded, lower=0, upper=1)],
"laplace": True, "laplace": True,
"Y": self.binary_Y, "Y": self.binary_Y,
"ep": True "ep": True
}, },
"Exponential_default": { #"Exponential_default": {
"model": GPy.likelihoods.exponential(), #"model": GPy.likelihoods.exponential(),
"link_f_constraints": [constrain_positive], #"link_f_constraints": [constrain_positive],
"Y": self.positive_Y, #"Y": self.positive_Y,
"laplace": True, #"laplace": True,
}, #},
"Poisson_default": { #"Poisson_default": {
"model": GPy.likelihoods.poisson(), #"model": GPy.likelihoods.poisson(),
"link_f_constraints": [constrain_positive], #"link_f_constraints": [constrain_positive],
"Y": self.integer_Y, #"Y": self.integer_Y,
"laplace": True, #"laplace": True,
"ep": False #Should work though... #"ep": False #Should work though...
}, #},
"Gamma_default": { #"Gamma_default": {
"model": GPy.likelihoods.gamma(), #"model": GPy.likelihoods.gamma(),
"link_f_constraints": [constrain_positive], #"link_f_constraints": [constrain_positive],
"Y": self.positive_Y, #"Y": self.positive_Y,
"laplace": True #"laplace": True
} #}
} }
for name, attributes in noise_models.iteritems(): for name, attributes in noise_models.iteritems():
@ -286,8 +300,8 @@ class TestNoiseModels(object):
else: else:
ep = False ep = False
if len(param_vals) > 1: #if len(param_vals) > 1:
raise NotImplementedError("Cannot support multiple params in likelihood yet!") #raise NotImplementedError("Cannot support multiple params in likelihood yet!")
#Required by all #Required by all
#Normal derivatives #Normal derivatives
@ -302,13 +316,13 @@ class TestNoiseModels(object):
yield self.t_d3logpdf_df3, model, Y, f yield self.t_d3logpdf_df3, model, Y, f
yield self.t_d3logpdf_dlink3, model, Y, f, link_f_constraints yield self.t_d3logpdf_dlink3, model, Y, f, link_f_constraints
#Params #Params
yield self.t_dlogpdf_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_dlogpdf_df_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_df_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_d2logpdf2_df2_dparams, model, Y, f, param_vals, param_constraints yield self.t_d2logpdf2_df2_dparams, model, Y, f, param_vals, param_names, param_constraints
#Link params #Link params
yield self.t_dlogpdf_link_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_link_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_dlogpdf_dlink_dparams, model, Y, f, param_vals, param_constraints yield self.t_dlogpdf_dlink_dparams, model, Y, f, param_vals, param_names, param_constraints
yield self.t_d2logpdf2_dlink2_dparams, model, Y, f, param_vals, param_constraints yield self.t_d2logpdf2_dlink2_dparams, model, Y, f, param_vals, param_names, param_constraints
#laplace likelihood gradcheck #laplace likelihood gradcheck
yield self.t_laplace_fit_rbf_white, model, self.X, Y, f, self.step, param_vals, param_names, param_constraints yield self.t_laplace_fit_rbf_white, model, self.X, Y, f, self.step, param_vals, param_names, param_constraints
@ -370,33 +384,33 @@ class TestNoiseModels(object):
# df_dparams # # df_dparams #
############## ##############
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_dparams(self, model, Y, f, params, params_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.logpdf, model.dlogpdf_dtheta, dparam_checkgrad(model.logpdf, model.dlogpdf_dtheta,
params, args=(f, Y), constraints=param_constraints, params, params_names, args=(f, Y), constraints=param_constraints,
randomize=True, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_df_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_df_dparams(self, model, Y, f, params, params_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.dlogpdf_df, model.dlogpdf_df_dtheta, dparam_checkgrad(model.dlogpdf_df, model.dlogpdf_df_dtheta,
params, args=(f, Y), constraints=param_constraints, params, params_names, args=(f, Y), constraints=param_constraints,
randomize=True, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_d2logpdf2_df2_dparams(self, model, Y, f, params, param_constraints): def t_d2logpdf2_df2_dparams(self, model, Y, f, params, params_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.d2logpdf_df2, model.d2logpdf_df2_dtheta, dparam_checkgrad(model.d2logpdf_df2, model.d2logpdf_df2_dtheta,
params, args=(f, Y), constraints=param_constraints, params, params_names, args=(f, Y), constraints=param_constraints,
randomize=True, verbose=True) randomize=False, verbose=True)
) )
################ ################
@ -454,32 +468,32 @@ class TestNoiseModels(object):
# dlink_dparams # # dlink_dparams #
################# #################
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_link_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_link_dparams(self, model, Y, f, params, param_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.logpdf_link, model.dlogpdf_link_dtheta, dparam_checkgrad(model.logpdf_link, model.dlogpdf_link_dtheta,
params, args=(f, Y), constraints=param_constraints, params, param_names, args=(f, Y), constraints=param_constraints,
randomize=False, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_dlogpdf_dlink_dparams(self, model, Y, f, params, param_constraints): def t_dlogpdf_dlink_dparams(self, model, Y, f, params, param_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.dlogpdf_dlink, model.dlogpdf_dlink_dtheta, dparam_checkgrad(model.dlogpdf_dlink, model.dlogpdf_dlink_dtheta,
params, args=(f, Y), constraints=param_constraints, params, param_names, args=(f, Y), constraints=param_constraints,
randomize=False, verbose=True) randomize=False, verbose=True)
) )
@with_setup(setUp, tearDown) @with_setup(setUp, tearDown)
def t_d2logpdf2_dlink2_dparams(self, model, Y, f, params, param_constraints): def t_d2logpdf2_dlink2_dparams(self, model, Y, f, params, param_names, param_constraints):
print "\n{}".format(inspect.stack()[0][3]) print "\n{}".format(inspect.stack()[0][3])
print model print model
assert ( assert (
dparam_checkgrad(model.d2logpdf_dlink2, model.d2logpdf_dlink2_dtheta, dparam_checkgrad(model.d2logpdf_dlink2, model.d2logpdf_dlink2_dtheta,
params, args=(f, Y), constraints=param_constraints, params, param_names, args=(f, Y), constraints=param_constraints,
randomize=False, verbose=True) randomize=False, verbose=True)
) )
@ -493,18 +507,23 @@ class TestNoiseModels(object):
Y = Y/Y.max() Y = Y/Y.max()
white_var = 1e-6 white_var = 1e-6
kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), model) laplace_likelihood = GPy.inference.latent_function_inference.LaplaceInference()
m = GPy.models.GPRegression(X.copy(), Y.copy(), kernel, likelihood=laplace_likelihood) m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=model, inference_method=laplace_likelihood)
m.ensure_default_constraints() m.ensure_default_constraints()
m.constrain_fixed('white', white_var) m['white'].constrain_fixed(white_var)
for param_num in range(len(param_names)): #Set constraints
name = param_names[param_num] for constrain_param, constraint in constraints:
m[name] = param_vals[param_num] constraint(constrain_param, m)
constraints[param_num](name, m)
print m print m
m.randomize() m.randomize()
#Set params
for param_num in range(len(param_names)):
name = param_names[param_num]
m[name] = param_vals[param_num]
#m.optimize(max_iters=8) #m.optimize(max_iters=8)
print m print m
m.checkgrad(verbose=1, step=step) m.checkgrad(verbose=1, step=step)
@ -525,10 +544,10 @@ class TestNoiseModels(object):
Y = Y/Y.max() Y = Y/Y.max()
white_var = 1e-6 white_var = 1e-6
kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
ep_likelihood = GPy.likelihoods.EP(Y.copy(), model) ep_inf = GPy.inference.latent_function_inference.EP()
m = GPy.models.GPRegression(X.copy(), Y.copy(), kernel, likelihood=ep_likelihood) m = GPy.core.GP(X.copy(), Y.copy(), kernel=kernel, likelihood=model, inference_method=ep_inf)
m.ensure_default_constraints() m.ensure_default_constraints()
m.constrain_fixed('white', white_var) m['white'].constrain_fixed(white_var)
for param_num in range(len(param_names)): for param_num in range(len(param_names)):
name = param_names[param_num] name = param_names[param_num]
@ -559,8 +578,8 @@ class LaplaceTests(unittest.TestCase):
self.var = 0.2 self.var = 0.2
self.var = np.random.rand(1) self.var = np.random.rand(1)
self.stu_t = GPy.likelihoods.student_t(deg_free=5, sigma2=self.var) self.stu_t = GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var)
self.gauss = GPy.likelihoods.gaussian(gp_transformations.Log(), variance=self.var, D=self.D, N=self.N) self.gauss = GPy.likelihoods.Gaussian(gp_link=link_functions.Log(), variance=self.var)
#Make a bigger step as lower bound can be quite curved #Make a bigger step as lower bound can be quite curved
self.step = 1e-6 self.step = 1e-6
@ -584,7 +603,7 @@ class LaplaceTests(unittest.TestCase):
noise = np.random.randn(*self.X.shape)*self.real_std noise = np.random.randn(*self.X.shape)*self.real_std
self.Y = np.sin(self.X*2*np.pi) + noise self.Y = np.sin(self.X*2*np.pi) + noise
self.f = np.random.rand(self.N, 1) self.f = np.random.rand(self.N, 1)
self.gauss = GPy.likelihoods.gaussian(variance=self.var, D=self.D, N=self.N) self.gauss = GPy.likelihoods.Gaussian(variance=self.var)
dlogpdf_df = functools.partial(self.gauss.dlogpdf_df, y=self.Y) dlogpdf_df = functools.partial(self.gauss.dlogpdf_df, y=self.Y)
d2logpdf_df2 = functools.partial(self.gauss.d2logpdf_df2, y=self.Y) d2logpdf_df2 = functools.partial(self.gauss.d2logpdf_df2, y=self.Y)
@ -605,23 +624,27 @@ class LaplaceTests(unittest.TestCase):
#Yc = Y.copy() #Yc = Y.copy()
#Yc[75:80] += 1 #Yc[75:80] += 1
kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
kernel2 = kernel1.copy() #FIXME: Make sure you can copy kernels when params is fixed
#kernel2 = kernel1.copy()
kernel2 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1) gauss_distr1 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
m1.constrain_fixed('white', 1e-6) exact_inf = GPy.inference.latent_function_inference.ExactGaussianInference()
m1['noise'] = initial_var_guess m1 = GPy.core.GP(X, Y.copy(), kernel=kernel1, likelihood=gauss_distr1, inference_method=exact_inf)
m1.constrain_bounded('noise', 1e-4, 10) m1['white'].constrain_fixed(1e-6)
m1.constrain_bounded('rbf', 1e-4, 10) m1['variance'] = initial_var_guess
m1['variance'].constrain_bounded(1e-4, 10)
m1['rbf'].constrain_bounded(1e-4, 10)
m1.ensure_default_constraints() m1.ensure_default_constraints()
m1.randomize() m1.randomize()
gauss_distr = GPy.likelihoods.gaussian(variance=initial_var_guess, D=1, N=Y.shape[0]) gauss_distr2 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), gauss_distr) laplace_inf = GPy.inference.latent_function_inference.LaplaceInference()
m2 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel2, likelihood=laplace_likelihood) m2 = GPy.core.GP(X, Y.copy(), kernel=kernel2, likelihood=gauss_distr2, inference_method=laplace_inf)
m2.ensure_default_constraints() m2.ensure_default_constraints()
m2.constrain_fixed('white', 1e-6) m2['white'].constrain_fixed(1e-6)
m2.constrain_bounded('rbf', 1e-4, 10) m2['rbf'].constrain_bounded(1e-4, 10)
m2.constrain_bounded('noise', 1e-4, 10) m2['variance'].constrain_bounded(1e-4, 10)
m2.randomize() m2.randomize()
if debug: if debug:
@ -667,7 +690,7 @@ class LaplaceTests(unittest.TestCase):
#Check Y's are the same #Check Y's are the same
np.testing.assert_almost_equal(Y, m2.likelihood.Y, decimal=5) np.testing.assert_almost_equal(m1.Y, m2.Y, decimal=5)
#Check marginals are the same #Check marginals are the same
np.testing.assert_almost_equal(m1.log_likelihood(), m2.log_likelihood(), decimal=2) np.testing.assert_almost_equal(m1.log_likelihood(), m2.log_likelihood(), decimal=2)
#Check marginals are the same with random #Check marginals are the same with random