mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-18 13:55:14 +02:00
[huge merge] the second
This commit is contained in:
Max Zwiessele
parent
180650ec85
commit
187f85c239
35 changed files with 40 additions and 3018 deletions
|
|
@ -1,406 +0,0 @@
|
|||
# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
#
|
||||
#Parts of this file were influenced by the Matlab GPML framework written by
|
||||
#Carl Edward Rasmussen & Hannes Nickisch, however all bugs are our own.
|
||||
#
|
||||
#The GPML code is released under the FreeBSD License.
|
||||
#Copyright (c) 2005-2013 Carl Edward Rasmussen & Hannes Nickisch. All rights reserved.
|
||||
#
|
||||
#The code and associated documentation is available from
|
||||
#http://gaussianprocess.org/gpml/code.
|
||||
|
||||
import numpy as np
|
||||
import scipy as sp
|
||||
from likelihood import likelihood
|
||||
from ..util.linalg import mdot, jitchol, pddet, dpotrs
|
||||
from functools import partial as partial_func
|
||||
import warnings
|
||||
|
||||
class Laplace(likelihood):
|
||||
"""Laplace approximation to a posterior"""
|
||||
|
||||
def __init__(self, data, noise_model, extra_data=None):
|
||||
"""
|
||||
Laplace Approximation
|
||||
|
||||
Find the moments \hat{f} and the hessian at this point
|
||||
(using Newton-Raphson) of the unnormalised posterior
|
||||
|
||||
Compute the GP variables (i.e. generate some Y^{squiggle} and
|
||||
z^{squiggle} which makes a gaussian the same as the laplace
|
||||
approximation to the posterior, but normalised
|
||||
|
||||
Arguments
|
||||
---------
|
||||
|
||||
:param data: array of data the likelihood function is approximating
|
||||
:type data: NxD
|
||||
:param noise_model: likelihood function - subclass of noise_model
|
||||
:type noise_model: noise_model
|
||||
:param extra_data: additional data used by some likelihood functions,
|
||||
"""
|
||||
self.data = data
|
||||
self.noise_model = noise_model
|
||||
self.extra_data = extra_data
|
||||
|
||||
#Inital values
|
||||
self.N, self.D = self.data.shape
|
||||
self.is_heteroscedastic = True
|
||||
self.Nparams = 0
|
||||
self.NORMAL_CONST = ((0.5 * self.N) * np.log(2 * np.pi))
|
||||
|
||||
self.restart()
|
||||
likelihood.__init__(self)
|
||||
|
||||
def restart(self):
|
||||
"""
|
||||
Reset likelihood variables to their defaults
|
||||
"""
|
||||
#Initial values for the GP variables
|
||||
self.Y = np.zeros((self.N, 1))
|
||||
self.covariance_matrix = np.eye(self.N)
|
||||
self.precision = np.ones(self.N)[:, None]
|
||||
self.Z = 0
|
||||
self.YYT = None
|
||||
|
||||
self.old_Ki_f = None
|
||||
self.bad_fhat = False
|
||||
|
||||
def predictive_values(self,mu,var,full_cov,**noise_args):
|
||||
if full_cov:
|
||||
raise NotImplementedError, "Cannot make correlated predictions with an EP likelihood"
|
||||
return self.noise_model.predictive_values(mu,var,**noise_args)
|
||||
|
||||
def log_predictive_density(self, y_test, mu_star, var_star):
|
||||
"""
|
||||
Calculation of the log predictive density
|
||||
|
||||
.. math:
|
||||
p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*})
|
||||
|
||||
:param y_test: test observations (y_{*})
|
||||
:type y_test: (Nx1) array
|
||||
:param mu_star: predictive mean of gaussian p(f_{*}|mu_{*}, var_{*})
|
||||
:type mu_star: (Nx1) array
|
||||
:param var_star: predictive variance of gaussian p(f_{*}|mu_{*}, var_{*})
|
||||
:type var_star: (Nx1) array
|
||||
"""
|
||||
return self.noise_model.log_predictive_density(y_test, mu_star, var_star)
|
||||
|
||||
def _get_params(self):
|
||||
return np.asarray(self.noise_model._get_params())
|
||||
|
||||
def _get_param_names(self):
|
||||
return self.noise_model._get_param_names()
|
||||
|
||||
def _set_params(self, p):
|
||||
return self.noise_model._set_params(p)
|
||||
|
||||
def _shared_gradients_components(self):
|
||||
d3lik_d3fhat = self.noise_model.d3logpdf_df3(self.f_hat, self.data, extra_data=self.extra_data)
|
||||
dL_dfhat = 0.5*(np.diag(self.Ki_W_i)[:, None]*d3lik_d3fhat).T #why isn't this -0.5?
|
||||
I_KW_i = np.eye(self.N) - np.dot(self.K, self.Wi_K_i)
|
||||
return dL_dfhat, I_KW_i
|
||||
|
||||
def _Kgradients(self):
|
||||
"""
|
||||
Gradients with respect to prior kernel parameters dL_dK to be chained
|
||||
with dK_dthetaK to give dL_dthetaK
|
||||
:returns: dL_dK matrix
|
||||
:rtype: Matrix (1 x num_kernel_params)
|
||||
"""
|
||||
dL_dfhat, I_KW_i = self._shared_gradients_components()
|
||||
dlp = self.noise_model.dlogpdf_df(self.f_hat, self.data, extra_data=self.extra_data)
|
||||
|
||||
#Explicit
|
||||
#expl_a = np.dot(self.Ki_f, self.Ki_f.T)
|
||||
#expl_b = self.Wi_K_i
|
||||
#expl = 0.5*expl_a - 0.5*expl_b
|
||||
#dL_dthetaK_exp = dK_dthetaK(expl, X)
|
||||
|
||||
#Implicit
|
||||
impl = mdot(dlp, dL_dfhat, I_KW_i)
|
||||
|
||||
#No longer required as we are computing these in the gp already
|
||||
#otherwise we would take them away and add them back
|
||||
#dL_dthetaK_imp = dK_dthetaK(impl, X)
|
||||
#dL_dthetaK = dL_dthetaK_exp + dL_dthetaK_imp
|
||||
#dL_dK = expl + impl
|
||||
|
||||
#No need to compute explicit as we are computing dZ_dK to account
|
||||
#for the difference between the K gradients of a normal GP,
|
||||
#and the K gradients including the implicit part
|
||||
dL_dK = impl
|
||||
return dL_dK
|
||||
|
||||
def _gradients(self, partial):
|
||||
"""
|
||||
Gradients with respect to likelihood parameters (dL_dthetaL)
|
||||
|
||||
:param partial: Not needed by this likelihood
|
||||
:type partial: lambda function
|
||||
: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 = self.noise_model._laplace_gradients(self.f_hat, self.data, extra_data=self.extra_data)
|
||||
|
||||
#len(dlik_dthetaL)
|
||||
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
|
||||
dL_dthetaL_exp = ( np.sum(dlik_dthetaL[:, thetaL_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
|
||||
dfhat_dthetaL = mdot(I_KW_i, self.K, dlik_grad_dthetaL[:, thetaL_i])
|
||||
dL_dthetaL_imp = np.dot(dL_dfhat, dfhat_dthetaL)
|
||||
dL_dthetaL[thetaL_i] = dL_dthetaL_exp + dL_dthetaL_imp
|
||||
|
||||
return dL_dthetaL
|
||||
|
||||
def _compute_GP_variables(self):
|
||||
"""
|
||||
Generate data Y which would give the normal distribution identical
|
||||
to the laplace approximation to the posterior, but normalised
|
||||
|
||||
GPy expects a likelihood to be gaussian, so need to caluclate
|
||||
the data Y^{\tilde} that makes the posterior match that found
|
||||
by a laplace approximation to a non-gaussian likelihood but with
|
||||
a gaussian likelihood
|
||||
|
||||
Firstly,
|
||||
The hessian of the unormalised posterior distribution is (K^{-1} + W)^{-1},
|
||||
i.e. z*N(f|f^{\hat}, (K^{-1} + W)^{-1}) but this assumes a non-gaussian likelihood,
|
||||
we wish to find the hessian \Sigma^{\tilde}
|
||||
that has the same curvature but using our new simulated data Y^{\tilde}
|
||||
i.e. we do N(Y^{\tilde}|f^{\hat}, \Sigma^{\tilde})N(f|0, K) = z*N(f|f^{\hat}, (K^{-1} + W)^{-1})
|
||||
and we wish to find what Y^{\tilde} and \Sigma^{\tilde}
|
||||
We find that Y^{\tilde} = W^{-1}(K^{-1} + W)f^{\hat} and \Sigma^{tilde} = W^{-1}
|
||||
|
||||
Secondly,
|
||||
GPy optimizes the log marginal log p(y) = -0.5*ln|K+\Sigma^{\tilde}| - 0.5*Y^{\tilde}^{T}(K^{-1} + \Sigma^{tilde})^{-1}Y + lik.Z
|
||||
So we can suck up any differences between that and our log marginal likelihood approximation
|
||||
p^{\squiggle}(y) = -0.5*f^{\hat}K^{-1}f^{\hat} + log p(y|f^{\hat}) - 0.5*log |K||K^{-1} + W|
|
||||
which we want to optimize instead, by equating them and rearranging, the difference is added onto
|
||||
the log p(y) that GPy optimizes by default
|
||||
|
||||
Thirdly,
|
||||
Since we have gradients that depend on how we move f^{\hat}, we have implicit components
|
||||
aswell as the explicit dL_dK, we hold these differences in dZ_dK and add them to dL_dK in the
|
||||
gp.py code
|
||||
"""
|
||||
Wi = 1.0/self.W
|
||||
self.Sigma_tilde = np.diagflat(Wi)
|
||||
|
||||
Y_tilde = Wi*self.Ki_f + self.f_hat
|
||||
|
||||
self.Wi_K_i = self.W12BiW12
|
||||
ln_det_Wi_K = pddet(self.Sigma_tilde + self.K)
|
||||
lik = self.noise_model.logpdf(self.f_hat, self.data, extra_data=self.extra_data)
|
||||
y_Wi_K_i_y = mdot(Y_tilde.T, self.Wi_K_i, Y_tilde)
|
||||
|
||||
Z_tilde = (+ lik
|
||||
- 0.5*self.ln_B_det
|
||||
+ 0.5*ln_det_Wi_K
|
||||
- 0.5*self.f_Ki_f
|
||||
+ 0.5*y_Wi_K_i_y
|
||||
+ self.NORMAL_CONST
|
||||
)
|
||||
|
||||
#Convert to float as its (1, 1) and Z must be a scalar
|
||||
self.Z = np.float64(Z_tilde)
|
||||
self.Y = Y_tilde
|
||||
self.YYT = np.dot(self.Y, self.Y.T)
|
||||
self.covariance_matrix = self.Sigma_tilde
|
||||
self.precision = 1.0 / np.diag(self.covariance_matrix)[:, None]
|
||||
|
||||
#Compute dZ_dK which is how the approximated distributions gradients differ from the dL_dK computed for other likelihoods
|
||||
self.dZ_dK = self._Kgradients()
|
||||
#+ 0.5*self.Wi_K_i - 0.5*np.dot(self.Ki_f, self.Ki_f.T) #since we are not adding the K gradients explicit part theres no need to compute this again
|
||||
|
||||
def fit_full(self, K):
|
||||
"""
|
||||
The laplace approximation algorithm, find K and expand hessian
|
||||
For nomenclature see Rasmussen & Williams 2006 - modified for numerical stability
|
||||
|
||||
:param K: Prior covariance matrix evaluated at locations X
|
||||
:type K: NxN matrix
|
||||
"""
|
||||
self.K = K.copy()
|
||||
|
||||
#Find mode
|
||||
self.f_hat = self.rasm_mode(self.K)
|
||||
|
||||
#Compute hessian and other variables at mode
|
||||
self._compute_likelihood_variables()
|
||||
|
||||
#Compute fake variables replicating laplace approximation to posterior
|
||||
self._compute_GP_variables()
|
||||
|
||||
def _compute_likelihood_variables(self):
|
||||
"""
|
||||
Compute the variables required to compute gaussian Y variables
|
||||
"""
|
||||
#At this point get the hessian matrix (or vector as W is diagonal)
|
||||
self.W = -self.noise_model.d2logpdf_df2(self.f_hat, self.data, extra_data=self.extra_data)
|
||||
|
||||
if not self.noise_model.log_concave:
|
||||
i = self.W < 1e-6
|
||||
if np.any(i):
|
||||
warnings.warn('truncating non log-concave likelihood curvature')
|
||||
# FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
|
||||
self.W[i] = 1e-6
|
||||
|
||||
self.W12BiW12, self.ln_B_det = self._compute_B_statistics(self.K, self.W, np.eye(self.N))
|
||||
|
||||
self.Ki_f = self.Ki_f
|
||||
self.f_Ki_f = np.dot(self.f_hat.T, self.Ki_f)
|
||||
self.Ki_W_i = self.K - mdot(self.K, self.W12BiW12, self.K)
|
||||
|
||||
def _compute_B_statistics(self, K, W, a):
|
||||
"""
|
||||
Rasmussen suggests the use of a numerically stable positive definite matrix B
|
||||
Which has a positive diagonal element and can be easyily inverted
|
||||
|
||||
:param K: Prior Covariance matrix evaluated at locations X
|
||||
:type K: NxN matrix
|
||||
:param W: Negative hessian at a point (diagonal matrix)
|
||||
:type W: Vector of diagonal values of hessian (1xN)
|
||||
:param a: Matrix to calculate W12BiW12a
|
||||
:type a: Matrix NxN
|
||||
:returns: (W12BiW12a, ln_B_det)
|
||||
"""
|
||||
if not self.noise_model.log_concave:
|
||||
#print "Under 1e-10: {}".format(np.sum(W < 1e-6))
|
||||
W[W < 1e-10] = 1e-10 # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
|
||||
# If the likelihood is non-log-concave. We wan't to say that there is a negative variance
|
||||
# To cause the posterior to become less certain than the prior and likelihood,
|
||||
# This is a property only held by non-log-concave likelihoods
|
||||
|
||||
|
||||
#W is diagonal so its sqrt is just the sqrt of the diagonal elements
|
||||
W_12 = np.sqrt(W)
|
||||
B = np.eye(self.N) + W_12*K*W_12.T
|
||||
L = jitchol(B)
|
||||
|
||||
W12BiW12a = W_12*dpotrs(L, np.asfortranarray(W_12*a), lower=1)[0]
|
||||
ln_B_det = 2*np.sum(np.log(np.diag(L)))
|
||||
return W12BiW12a, ln_B_det
|
||||
|
||||
def rasm_mode(self, K, MAX_ITER=40):
|
||||
"""
|
||||
Rasmussen's numerically stable mode finding
|
||||
For nomenclature see Rasmussen & Williams 2006
|
||||
Influenced by GPML (BSD) code, all errors are our own
|
||||
|
||||
:param K: Covariance matrix evaluated at locations X
|
||||
:type K: NxD matrix
|
||||
:param MAX_ITER: Maximum number of iterations of newton-raphson before forcing finish of optimisation
|
||||
:type MAX_ITER: scalar
|
||||
:returns: f_hat, mode on which to make laplace approxmiation
|
||||
:rtype: NxD matrix
|
||||
"""
|
||||
#old_Ki_f = np.zeros((self.N, 1))
|
||||
|
||||
#Start f's at zero originally of if we have gone off track, try restarting
|
||||
if self.old_Ki_f is None or self.bad_fhat:
|
||||
old_Ki_f = np.random.rand(self.N, 1)/50.0
|
||||
#old_Ki_f = self.Y
|
||||
f = np.dot(K, old_Ki_f)
|
||||
else:
|
||||
#Start at the old best point
|
||||
old_Ki_f = self.old_Ki_f.copy()
|
||||
f = self.f_hat.copy()
|
||||
|
||||
new_obj = -np.inf
|
||||
old_obj = np.inf
|
||||
|
||||
def obj(Ki_f, f):
|
||||
return -0.5*np.dot(Ki_f.T, f) + self.noise_model.logpdf(f, self.data, extra_data=self.extra_data)
|
||||
|
||||
difference = np.inf
|
||||
epsilon = 1e-7
|
||||
#step_size = 1
|
||||
#rs = 0
|
||||
i = 0
|
||||
|
||||
while difference > epsilon and i < MAX_ITER:
|
||||
W = -self.noise_model.d2logpdf_df2(f, self.data, extra_data=self.extra_data)
|
||||
|
||||
W_f = W*f
|
||||
grad = self.noise_model.dlogpdf_df(f, self.data, extra_data=self.extra_data)
|
||||
|
||||
b = W_f + grad
|
||||
W12BiW12Kb, _ = self._compute_B_statistics(K, W.copy(), np.dot(K, b))
|
||||
|
||||
#Work out the DIRECTION that we want to move in, but don't choose the stepsize yet
|
||||
full_step_Ki_f = b - W12BiW12Kb
|
||||
dKi_f = full_step_Ki_f - old_Ki_f
|
||||
|
||||
f_old = f.copy()
|
||||
def inner_obj(step_size, old_Ki_f, dKi_f, K):
|
||||
Ki_f = old_Ki_f + step_size*dKi_f
|
||||
f = np.dot(K, Ki_f)
|
||||
# This is nasty, need to set something within an optimization though
|
||||
self.tmp_Ki_f = Ki_f.copy()
|
||||
self.tmp_f = f.copy()
|
||||
return -obj(Ki_f, f)
|
||||
|
||||
i_o = partial_func(inner_obj, old_Ki_f=old_Ki_f, dKi_f=dKi_f, K=K)
|
||||
#Find the stepsize that minimizes the objective function using a brent line search
|
||||
#The tolerance and maxiter matter for speed! Seems to be best to keep them low and make more full
|
||||
#steps than get this exact then make a step, if B was bigger it might be the other way around though
|
||||
#new_obj = sp.optimize.minimize_scalar(i_o, method='brent', tol=1e-4, options={'maxiter':5}).fun
|
||||
new_obj = sp.optimize.brent(i_o, tol=1e-4, maxiter=10)
|
||||
f = self.tmp_f.copy()
|
||||
Ki_f = self.tmp_Ki_f.copy()
|
||||
|
||||
#Optimize without linesearch
|
||||
#f_old = f.copy()
|
||||
#update_passed = False
|
||||
#while not update_passed:
|
||||
#Ki_f = old_Ki_f + step_size*dKi_f
|
||||
#f = np.dot(K, Ki_f)
|
||||
|
||||
#old_obj = new_obj
|
||||
#new_obj = obj(Ki_f, f)
|
||||
#difference = new_obj - old_obj
|
||||
##print "difference: ",difference
|
||||
#if difference < 0:
|
||||
##print "Objective function rose", np.float(difference)
|
||||
##If the objective function isn't rising, restart optimization
|
||||
#step_size *= 0.8
|
||||
##print "Reducing step-size to {ss:.3} and restarting optimization".format(ss=step_size)
|
||||
##objective function isn't increasing, try reducing step size
|
||||
#f = f_old.copy() #it's actually faster not to go back to old location and just zigzag across the mode
|
||||
#old_obj = new_obj
|
||||
#rs += 1
|
||||
#else:
|
||||
#update_passed = True
|
||||
|
||||
#old_Ki_f = self.Ki_f.copy()
|
||||
|
||||
#difference = abs(new_obj - old_obj)
|
||||
#old_obj = new_obj.copy()
|
||||
difference = np.abs(np.sum(f - f_old)) + np.abs(np.sum(Ki_f - old_Ki_f))
|
||||
#difference = np.abs(np.sum(Ki_f - old_Ki_f))/np.float(self.N)
|
||||
old_Ki_f = Ki_f.copy()
|
||||
i += 1
|
||||
|
||||
self.old_Ki_f = old_Ki_f.copy()
|
||||
|
||||
#Warn of bad fits
|
||||
if difference > epsilon:
|
||||
self.bad_fhat = True
|
||||
warnings.warn("Not perfect f_hat fit difference: {}".format(difference))
|
||||
elif self.bad_fhat:
|
||||
self.bad_fhat = False
|
||||
warnings.warn("f_hat now perfect again")
|
||||
|
||||
self.Ki_f = Ki_f
|
||||
return f
|
||||
|
|
@ -1,222 +0,0 @@
|
|||
# Copyright (c) 2012, 2013 Ricardo Andrade
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
|
||||
import numpy as np
|
||||
from scipy import stats,special
|
||||
import scipy as sp
|
||||
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
|
||||
import gp_transformations
|
||||
from noise_distributions import NoiseDistribution
|
||||
|
||||
class Bernoulli(NoiseDistribution):
|
||||
"""
|
||||
Bernoulli likelihood
|
||||
|
||||
.. math::
|
||||
p(y_{i}|\\lambda(f_{i})) = \\lambda(f_{i})^{y_{i}}(1-f_{i})^{1-y_{i}}
|
||||
|
||||
.. Note::
|
||||
Y is expected to take values in {-1,1}
|
||||
Probit likelihood usually used
|
||||
"""
|
||||
def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False):
|
||||
super(Bernoulli, self).__init__(gp_link,analytical_mean,analytical_variance)
|
||||
if isinstance(gp_link , (gp_transformations.Heaviside, gp_transformations.Probit)):
|
||||
self.log_concave = True
|
||||
|
||||
def _preprocess_values(self,Y):
|
||||
"""
|
||||
Check if the values of the observations correspond to the values
|
||||
assumed by the likelihood function.
|
||||
|
||||
..Note:: Binary classification algorithm works better with classes {-1,1}
|
||||
"""
|
||||
Y_prep = Y.copy()
|
||||
Y1 = Y[Y.flatten()==1].size
|
||||
Y2 = Y[Y.flatten()==0].size
|
||||
assert Y1 + Y2 == Y.size, 'Bernoulli likelihood is meant to be used only with outputs in {0,1}.'
|
||||
Y_prep[Y.flatten() == 0] = -1
|
||||
return Y_prep
|
||||
|
||||
def _moments_match_analytical(self,data_i,tau_i,v_i):
|
||||
"""
|
||||
Moments match of the marginal approximation in EP algorithm
|
||||
|
||||
:param i: number of observation (int)
|
||||
:param tau_i: precision of the cavity distribution (float)
|
||||
:param v_i: mean/variance of the cavity distribution (float)
|
||||
"""
|
||||
if data_i == 1:
|
||||
sign = 1.
|
||||
elif data_i == 0:
|
||||
sign = -1
|
||||
else:
|
||||
raise ValueError("bad value for Bernouilli observation (0,1)")
|
||||
if isinstance(self.gp_link,gp_transformations.Probit):
|
||||
z = sign*v_i/np.sqrt(tau_i**2 + tau_i)
|
||||
Z_hat = std_norm_cdf(z)
|
||||
phi = std_norm_pdf(z)
|
||||
mu_hat = v_i/tau_i + sign*phi/(Z_hat*np.sqrt(tau_i**2 + tau_i))
|
||||
sigma2_hat = 1./tau_i - (phi/((tau_i**2+tau_i)*Z_hat))*(z+phi/Z_hat)
|
||||
|
||||
elif isinstance(self.gp_link,gp_transformations.Heaviside):
|
||||
a = sign*v_i/np.sqrt(tau_i)
|
||||
Z_hat = std_norm_cdf(a)
|
||||
N = std_norm_pdf(a)
|
||||
mu_hat = v_i/tau_i + sign*N/Z_hat/np.sqrt(tau_i)
|
||||
sigma2_hat = (1. - a*N/Z_hat - np.square(N/Z_hat))/tau_i
|
||||
if np.any(np.isnan([Z_hat, mu_hat, sigma2_hat])):
|
||||
stop
|
||||
else:
|
||||
raise ValueError("Exact moment matching not available for link {}".format(self.gp_link.gp_transformations.__name__))
|
||||
|
||||
return Z_hat, mu_hat, sigma2_hat
|
||||
|
||||
def _predictive_mean_analytical(self,mu,variance):
|
||||
|
||||
if isinstance(self.gp_link,gp_transformations.Probit):
|
||||
return stats.norm.cdf(mu/np.sqrt(1+variance))
|
||||
|
||||
elif isinstance(self.gp_link,gp_transformations.Heaviside):
|
||||
return stats.norm.cdf(mu/np.sqrt(variance))
|
||||
|
||||
else:
|
||||
raise NotImplementedError
|
||||
|
||||
def _predictive_variance_analytical(self,mu,variance, pred_mean):
|
||||
|
||||
if isinstance(self.gp_link,gp_transformations.Heaviside):
|
||||
return 0.
|
||||
else:
|
||||
raise NotImplementedError
|
||||
|
||||
def pdf_link(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Likelihood function given link(f)
|
||||
|
||||
.. math::
|
||||
p(y_{i}|\\lambda(f_{i})) = \\lambda(f_{i})^{y_{i}}(1-f_{i})^{1-y_{i}}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data not used in bernoulli
|
||||
:returns: likelihood evaluated for this point
|
||||
:rtype: float
|
||||
|
||||
.. Note:
|
||||
Each y_i must be in {0,1}
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
objective = (link_f**y) * ((1.-link_f)**(1.-y))
|
||||
return np.exp(np.sum(np.log(objective)))
|
||||
|
||||
def logpdf_link(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Log Likelihood function given link(f)
|
||||
|
||||
.. math::
|
||||
\\ln p(y_{i}|\\lambda(f_{i})) = y_{i}\\log\\lambda(f_{i}) + (1-y_{i})\\log (1-f_{i})
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data not used in bernoulli
|
||||
:returns: log likelihood evaluated at points link(f)
|
||||
:rtype: float
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
#objective = y*np.log(link_f) + (1.-y)*np.log(link_f)
|
||||
objective = np.where(y==1, np.log(link_f), np.log(1-link_f))
|
||||
return np.sum(objective)
|
||||
|
||||
def dlogpdf_dlink(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Gradient of the pdf at y, given link(f) w.r.t link(f)
|
||||
|
||||
.. math::
|
||||
\\frac{d\\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)} = \\frac{y_{i}}{\\lambda(f_{i})} - \\frac{(1 - y_{i})}{(1 - \\lambda(f_{i}))}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data not used in bernoulli
|
||||
:returns: gradient of log likelihood evaluated at points link(f)
|
||||
:rtype: Nx1 array
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
grad = (y/link_f) - (1.-y)/(1-link_f)
|
||||
return grad
|
||||
|
||||
def d2logpdf_dlink2(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Hessian at y, given link_f, w.r.t link_f the hessian will be 0 unless i == j
|
||||
i.e. second derivative logpdf at y given link(f_i) link(f_j) w.r.t link(f_i) and link(f_j)
|
||||
|
||||
|
||||
.. math::
|
||||
\\frac{d^{2}\\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)^{2}} = \\frac{-y_{i}}{\\lambda(f)^{2}} - \\frac{(1-y_{i})}{(1-\\lambda(f))^{2}}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data not used in bernoulli
|
||||
:returns: Diagonal of log hessian matrix (second derivative of log likelihood evaluated at points link(f))
|
||||
:rtype: Nx1 array
|
||||
|
||||
.. Note::
|
||||
Will return diagonal of hessian, since every where else it is 0, as the likelihood factorizes over cases
|
||||
(the distribution for y_i depends only on link(f_i) not on link(f_(j!=i))
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
d2logpdf_dlink2 = -y/(link_f**2) - (1-y)/((1-link_f)**2)
|
||||
return d2logpdf_dlink2
|
||||
|
||||
def d3logpdf_dlink3(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Third order derivative log-likelihood function at y given link(f) w.r.t link(f)
|
||||
|
||||
.. math::
|
||||
\\frac{d^{3} \\ln p(y_{i}|\\lambda(f_{i}))}{d^{3}\\lambda(f)} = \\frac{2y_{i}}{\\lambda(f)^{3}} - \\frac{2(1-y_{i}}{(1-\\lambda(f))^{3}}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data not used in bernoulli
|
||||
:returns: third derivative of log likelihood evaluated at points link(f)
|
||||
:rtype: Nx1 array
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
d3logpdf_dlink3 = 2*(y/(link_f**3) - (1-y)/((1-link_f)**3))
|
||||
return d3logpdf_dlink3
|
||||
|
||||
def _mean(self,gp):
|
||||
"""
|
||||
Mass (or density) function
|
||||
"""
|
||||
return self.gp_link.transf(gp)
|
||||
|
||||
def _variance(self,gp):
|
||||
"""
|
||||
Mass (or density) function
|
||||
"""
|
||||
p = self.gp_link.transf(gp)
|
||||
return p*(1.-p)
|
||||
|
||||
def samples(self, gp):
|
||||
"""
|
||||
Returns a set of samples of observations based on a given value of the latent variable.
|
||||
|
||||
:param gp: latent variable
|
||||
"""
|
||||
orig_shape = gp.shape
|
||||
gp = gp.flatten()
|
||||
ns = np.ones_like(gp, dtype=int)
|
||||
Ysim = np.random.binomial(ns, self.gp_link.transf(gp))
|
||||
return Ysim.reshape(orig_shape)
|
||||
|
|
@ -1,277 +0,0 @@
|
|||
# Copyright (c) 2012, 2013 Ricardo Andrade
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import numpy as np
|
||||
from scipy import stats, special
|
||||
import scipy as sp
|
||||
import gp_transformations
|
||||
from noise_distributions import NoiseDistribution
|
||||
from scipy import stats, integrate
|
||||
from scipy.special import gammaln, gamma
|
||||
|
||||
class StudentT(NoiseDistribution):
|
||||
"""
|
||||
Student T likelihood
|
||||
|
||||
For nomanclature see Bayesian Data Analysis 2003 p576
|
||||
|
||||
.. math::
|
||||
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):
|
||||
self.v = deg_free
|
||||
self.sigma2 = sigma2
|
||||
|
||||
self._set_params(np.asarray(sigma2))
|
||||
super(StudentT, self).__init__(gp_link,analytical_mean,analytical_variance)
|
||||
self.log_concave = False
|
||||
|
||||
def _get_params(self):
|
||||
return np.asarray(self.sigma2)
|
||||
|
||||
def _get_param_names(self):
|
||||
return ["t_noise_std2"]
|
||||
|
||||
def _set_params(self, x):
|
||||
self.sigma2 = float(x)
|
||||
|
||||
@property
|
||||
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):
|
||||
"""
|
||||
Likelihood function given link(f)
|
||||
|
||||
.. math::
|
||||
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} - \\lambda(f_{i}))^{2}}{\\sigma^{2}}\\right)\\right)^{\\frac{-v+1}{2}}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: likelihood evaluated for this point
|
||||
:rtype: float
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
#Careful gamma(big_number) is infinity!
|
||||
objective = ((np.exp(gammaln((self.v + 1)*0.5) - gammaln(self.v * 0.5))
|
||||
/ (np.sqrt(self.v * np.pi * self.sigma2)))
|
||||
* ((1 + (1./float(self.v))*((e**2)/float(self.sigma2)))**(-0.5*(self.v + 1)))
|
||||
)
|
||||
return np.prod(objective)
|
||||
|
||||
def logpdf_link(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Log Likelihood Function given link(f)
|
||||
|
||||
.. math::
|
||||
\\ln p(y_{i}|\lambda(f_{i})) = \\ln \\Gamma\\left(\\frac{v+1}{2}\\right) - \\ln \\Gamma\\left(\\frac{v}{2}\\right) - \\ln \\sqrt{v \\pi\\sigma^{2}} - \\frac{v+1}{2}\\ln \\left(1 + \\frac{1}{v}\\left(\\frac{(y_{i} - \lambda(f_{i}))^{2}}{\\sigma^{2}}\\right)\\right)
|
||||
|
||||
:param link_f: latent variables (link(f))
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: likelihood evaluated for this point
|
||||
:rtype: float
|
||||
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
objective = (+ gammaln((self.v + 1) * 0.5)
|
||||
- gammaln(self.v * 0.5)
|
||||
- 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))
|
||||
)
|
||||
return np.sum(objective)
|
||||
|
||||
def dlogpdf_dlink(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Gradient of the log likelihood function at y, given link(f) w.r.t link(f)
|
||||
|
||||
.. math::
|
||||
\\frac{d \\ln p(y_{i}|\lambda(f_{i}))}{d\\lambda(f)} = \\frac{(v+1)(y_{i}-\lambda(f_{i}))}{(y_{i}-\lambda(f_{i}))^{2} + \\sigma^{2}v}
|
||||
|
||||
:param link_f: latent variables (f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: gradient of likelihood evaluated at points
|
||||
:rtype: Nx1 array
|
||||
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
grad = ((self.v + 1) * e) / (self.v * self.sigma2 + (e**2))
|
||||
return grad
|
||||
|
||||
def d2logpdf_dlink2(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Hessian at y, given link(f), w.r.t link(f)
|
||||
i.e. second derivative logpdf at y given link(f_i) and link(f_j) w.r.t link(f_i) and link(f_j)
|
||||
The hessian will be 0 unless i == j
|
||||
|
||||
.. math::
|
||||
\\frac{d^{2} \\ln p(y_{i}|\lambda(f_{i}))}{d^{2}\\lambda(f)} = \\frac{(v+1)((y_{i}-\lambda(f_{i}))^{2} - \\sigma^{2}v)}{((y_{i}-\lambda(f_{i}))^{2} + \\sigma^{2}v)^{2}}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: Diagonal of hessian matrix (second derivative of likelihood evaluated at points f)
|
||||
:rtype: Nx1 array
|
||||
|
||||
.. Note::
|
||||
Will return diagonal of hessian, since every where else it is 0, as the likelihood factorizes over cases
|
||||
(the distribution for y_i depends only on link(f_i) not on link(f_(j!=i))
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
hess = ((self.v + 1)*(e**2 - self.v*self.sigma2)) / ((self.sigma2*self.v + e**2)**2)
|
||||
return hess
|
||||
|
||||
def d3logpdf_dlink3(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Third order derivative log-likelihood function at y given link(f) w.r.t link(f)
|
||||
|
||||
.. math::
|
||||
\\frac{d^{3} \\ln p(y_{i}|\lambda(f_{i}))}{d^{3}\\lambda(f)} = \\frac{-2(v+1)((y_{i} - \lambda(f_{i}))^3 - 3(y_{i} - \lambda(f_{i})) \\sigma^{2} v))}{((y_{i} - \lambda(f_{i})) + \\sigma^{2} v)^3}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: third derivative of likelihood evaluated at points f
|
||||
:rtype: Nx1 array
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
d3lik_dlink3 = ( -(2*(self.v + 1)*(-e)*(e**2 - 3*self.v*self.sigma2)) /
|
||||
((e**2 + self.sigma2*self.v)**3)
|
||||
)
|
||||
return d3lik_dlink3
|
||||
|
||||
def dlogpdf_link_dvar(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Gradient of the log-likelihood function at y given f, w.r.t variance parameter (t_noise)
|
||||
|
||||
.. math::
|
||||
\\frac{d \\ln p(y_{i}|\lambda(f_{i}))}{d\\sigma^{2}} = \\frac{v((y_{i} - \lambda(f_{i}))^{2} - \\sigma^{2})}{2\\sigma^{2}(\\sigma^{2}v + (y_{i} - \lambda(f_{i}))^{2})}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: derivative of likelihood evaluated at points f w.r.t variance parameter
|
||||
:rtype: float
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
dlogpdf_dvar = self.v*(e**2 - self.sigma2)/(2*self.sigma2*(self.sigma2*self.v + e**2))
|
||||
return np.sum(dlogpdf_dvar)
|
||||
|
||||
def dlogpdf_dlink_dvar(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Derivative of the dlogpdf_dlink w.r.t variance parameter (t_noise)
|
||||
|
||||
.. math::
|
||||
\\frac{d}{d\\sigma^{2}}(\\frac{d \\ln p(y_{i}|\lambda(f_{i}))}{df}) = \\frac{-2\\sigma v(v + 1)(y_{i}-\lambda(f_{i}))}{(y_{i}-\lambda(f_{i}))^2 + \\sigma^2 v)^2}
|
||||
|
||||
:param link_f: latent variables link_f
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: derivative of likelihood evaluated at points f w.r.t variance parameter
|
||||
:rtype: Nx1 array
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
dlogpdf_dlink_dvar = (self.v*(self.v+1)*(-e))/((self.sigma2*self.v + e**2)**2)
|
||||
return dlogpdf_dlink_dvar
|
||||
|
||||
def d2logpdf_dlink2_dvar(self, link_f, y, extra_data=None):
|
||||
"""
|
||||
Gradient of the hessian (d2logpdf_dlink2) w.r.t variance parameter (t_noise)
|
||||
|
||||
.. math::
|
||||
\\frac{d}{d\\sigma^{2}}(\\frac{d^{2} \\ln p(y_{i}|\lambda(f_{i}))}{d^{2}f}) = \\frac{v(v+1)(\\sigma^{2}v - 3(y_{i} - \lambda(f_{i}))^{2})}{(\\sigma^{2}v + (y_{i} - \lambda(f_{i}))^{2})^{3}}
|
||||
|
||||
:param link_f: latent variables link(f)
|
||||
:type link_f: Nx1 array
|
||||
:param y: data
|
||||
:type y: Nx1 array
|
||||
:param extra_data: extra_data which is not used in student t distribution
|
||||
:returns: derivative of hessian evaluated at points f and f_j w.r.t variance parameter
|
||||
:rtype: Nx1 array
|
||||
"""
|
||||
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
|
||||
e = y - link_f
|
||||
d2logpdf_dlink2_dvar = ( (self.v*(self.v+1)*(self.sigma2*self.v - 3*(e**2)))
|
||||
/ ((self.sigma2*self.v + (e**2))**3)
|
||||
)
|
||||
return d2logpdf_dlink2_dvar
|
||||
|
||||
def dlogpdf_link_dtheta(self, f, y, extra_data=None):
|
||||
dlogpdf_dvar = self.dlogpdf_link_dvar(f, y, extra_data=extra_data)
|
||||
return np.asarray([[dlogpdf_dvar]])
|
||||
|
||||
def dlogpdf_dlink_dtheta(self, f, y, extra_data=None):
|
||||
dlogpdf_dlink_dvar = self.dlogpdf_dlink_dvar(f, y, extra_data=extra_data)
|
||||
return dlogpdf_dlink_dvar
|
||||
|
||||
def d2logpdf_dlink2_dtheta(self, f, y, extra_data=None):
|
||||
d2logpdf_dlink2_dvar = self.d2logpdf_dlink2_dvar(f, y, extra_data=extra_data)
|
||||
return d2logpdf_dlink2_dvar
|
||||
|
||||
def _predictive_variance_analytical(self, mu, sigma, predictive_mean=None):
|
||||
"""
|
||||
Compute predictive variance of student_t*normal p(y*|f*)p(f*)
|
||||
|
||||
Need to find what the variance is at the latent points for a student t*normal p(y*|f*)p(f*)
|
||||
(((g((v+1)/2))/(g(v/2)*s*sqrt(v*pi)))*(1+(1/v)*((y-f)/s)^2)^(-(v+1)/2))
|
||||
*((1/(s*sqrt(2*pi)))*exp(-(1/(2*(s^2)))*((y-f)^2)))
|
||||
"""
|
||||
|
||||
#FIXME: Not correct
|
||||
#We want the variance around test points y which comes from int p(y*|f*)p(f*) df*
|
||||
#Var(y*) = Var(E[y*|f*]) + E[Var(y*|f*)]
|
||||
#Since we are given f* (mu) which is our mean (expected) value of y*|f* then the variance is the variance around this
|
||||
#Which was also given to us as (var)
|
||||
#We also need to know the expected variance of y* around samples f*, this is the variance of the student t distribution
|
||||
#However the variance of the student t distribution is not dependent on f, only on sigma and the degrees of freedom
|
||||
true_var = 1/(1/sigma**2 + 1/self.variance)
|
||||
|
||||
return true_var
|
||||
|
||||
def _predictive_mean_analytical(self, mu, sigma):
|
||||
"""
|
||||
Compute mean of the prediction
|
||||
"""
|
||||
#FIXME: Not correct
|
||||
return mu
|
||||
|
||||
def samples(self, gp):
|
||||
"""
|
||||
Returns a set of samples of observations based on a given value of the latent variable.
|
||||
|
||||
:param gp: latent variable
|
||||
"""
|
||||
orig_shape = gp.shape
|
||||
gp = gp.flatten()
|
||||
#FIXME: Very slow as we are computing a new random variable per input!
|
||||
#Can't get it to sample all at the same time
|
||||
#student_t_samples = np.array([stats.t.rvs(self.v, self.gp_link.transf(gpj),scale=np.sqrt(self.sigma2), size=1) for gpj in gp])
|
||||
dfs = np.ones_like(gp)*self.v
|
||||
scales = np.ones_like(gp)*np.sqrt(self.sigma2)
|
||||
student_t_samples = stats.t.rvs(dfs, loc=self.gp_link.transf(gp),
|
||||
scale=scales)
|
||||
return student_t_samples.reshape(orig_shape)
|
||||
Loading…
Add table
Add a link
Reference in a new issue