tidying in likelihoods

This commit is contained in:
James Hensman 2014-03-10 08:34:03 +00:00
parent 2876e5a07a
commit 5c7caef885
4 changed files with 23 additions and 53 deletions

View file

@ -60,10 +60,9 @@ def student_t_approx(optimize=True, plot=True):
m3['.*white'].constrain_fixed(1e-5) m3['.*white'].constrain_fixed(1e-5)
m3.randomize() m3.randomize()
debug = True debug = True
print m3
if debug: #TODO: remove
m3.optimize(messages=1) return m3
return m3
#Student t GP model on corrupt data #Student t GP model on corrupt data
t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd) t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)

View file

@ -51,16 +51,16 @@ class Laplace(object):
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)
self.f_hat = f_hat self.f_hat = f_hat
#Compute hessian and other variables at mode #Compute hessian and other variables at mode
log_marginal, woodbury_vector, woodbury_inv, dL_dK, dL_dthetaL = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, kern, Y_metadata) log_marginal, woodbury_inv, dL_dK, dL_dthetaL = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, kern, Y_metadata)
kern.update_gradients_full(dL_dK, X) kern.update_gradients_full(dL_dK, X)
likelihood.update_gradients(dL_dthetaL) 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=woodbury_inv, K=K), log_marginal, {'dL_dK':dL_dK} return Posterior(woodbury_vector=Ki_fhat, 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):
""" """
@ -137,7 +137,6 @@ class Laplace(object):
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
woodbury_vector : variable required for calculating the approximation to the covariance matrix
woodbury_inv : 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_kernel_params)
dL_dthetaL : array of derivatives (1 x num_likelihood_params) dL_dthetaL : array of derivatives (1 x num_likelihood_params)
@ -156,7 +155,6 @@ class Laplace(object):
#Compute vival matrices for derivatives #Compute vival matrices for derivatives
dW_df = -likelihood.d3logpdf_df3(f_hat, Y, Y_metadata=Y_metadata) # -d3lik_d3fhat dW_df = -likelihood.d3logpdf_df3(f_hat, Y, Y_metadata=Y_metadata) # -d3lik_d3fhat
woodbury_vector = likelihood.dlogpdf_df(f_hat, Y, Y_metadata=Y_metadata)
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. 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) #BiK, _ = dpotrs(L, K, lower=1)
#dL_dfhat = 0.5*np.diag(BiK)[:, None]*dW_df #dL_dfhat = 0.5*np.diag(BiK)[:, None]*dW_df
@ -170,7 +168,7 @@ class Laplace(object):
explicit_part = 0.5*(np.dot(Ki_f, Ki_f.T) - K_Wi_i) explicit_part = 0.5*(np.dot(Ki_f, Ki_f.T) - K_Wi_i)
#Implicit #Implicit
implicit_part = np.dot(woodbury_vector, dL_dfhat.T).dot(I_KW_i) implicit_part = np.dot(Ki_f, dL_dfhat.T).dot(I_KW_i)
dL_dK = explicit_part + implicit_part dL_dK = explicit_part + implicit_part
else: else:
@ -201,7 +199,7 @@ class Laplace(object):
else: else:
dL_dthetaL = np.zeros(likelihood.size) dL_dthetaL = np.zeros(likelihood.size)
return log_marginal, woodbury_vector, K_Wi_i, dL_dK, dL_dthetaL return log_marginal, K_Wi_i, dL_dK, dL_dthetaL
def _compute_B_statistics(self, K, W, log_concave): def _compute_B_statistics(self, K, W, log_concave):
""" """

View file

@ -120,7 +120,7 @@ class Likelihood(Parameterized):
return z, mean, variance return z, mean, variance
def _predictive_mean(self,mu,variance): def _predictive_mean(self, mu, variance):
""" """
Quadrature calculation of the predictive mean: E(Y_star|Y) = E( E(Y_star|f_star, Y) ) Quadrature calculation of the predictive mean: E(Y_star|Y) = E( E(Y_star|f_star, Y) )
@ -135,6 +135,10 @@ class Likelihood(Parameterized):
return mean return mean
def _conditional_mean(self, f):
"""Quadrature calculation of the conditional mean: E(Y_star|f)"""
raise NotImplementedError, "implement this function to make predictions"
def _predictive_variance(self,mu,variance,predictive_mean=None): def _predictive_variance(self,mu,variance,predictive_mean=None):
""" """
Numerical approximation to the predictive variance: V(Y_star) Numerical approximation to the predictive variance: V(Y_star)
@ -358,50 +362,16 @@ class Likelihood(Parameterized):
return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta
def predictive_values(self, mu, var, full_cov=False, sampling=False, num_samples=10000): def predictive_values(self, mu, var):
""" """
Compute mean, variance and conficence interval (percentiles 5 and 95) of the prediction. Compute mean, variance of the predictive distibution.
:param mu: mean of the latent variable, f, of posterior :param mu: mean of the latent variable, f, of posterior
:param var: variance of the latent variable, f, of posterior :param var: variance of the latent variable, f, of posterior
:param full_cov: whether to use the full covariance or just the diagonal
:type full_cov: Boolean
:param num_samples: number of samples to use in computing quantiles and
possibly mean variance
:type num_samples: integer
:param sampling: Whether to use samples for mean and variances anyway
:type sampling: Boolean
""" """
pred_mean = self.predictive_mean(mu, var)
if sampling: pred_var = self.predictive_variance(mu, var, pred_mean)
#Get gp_samples f* using posterior mean and variance return pred_mean, pred_var
if not full_cov:
gp_samples = np.random.multivariate_normal(mu.flatten(), np.diag(var.flatten()),
size=num_samples).T
else:
gp_samples = np.random.multivariate_normal(mu.flatten(), var,
size=num_samples).T
#Push gp samples (f*) through likelihood to give p(y*|f*)
samples = self.samples(gp_samples)
axis=-1
#Calculate mean, variance and precentiles from samples
print "WARNING: Using sampling to calculate mean, variance and predictive quantiles."
pred_mean = np.mean(samples, axis=axis)[:,None]
pred_var = np.var(samples, axis=axis)[:,None]
q1 = np.percentile(samples, 2.5, axis=axis)[:,None]
q3 = np.percentile(samples, 97.5, axis=axis)[:,None]
else:
pred_mean = self.predictive_mean(mu, var)
pred_var = self.predictive_variance(mu, var, pred_mean)
print "WARNING: Predictive quantiles are only computed when sampling."
q1 = np.repeat(np.nan,pred_mean.size)[:,None]
q3 = q1.copy()
return pred_mean, pred_var, q1, q3
def samples(self, gp): def samples(self, gp):
""" """

View file

@ -25,8 +25,11 @@ class Poisson(Likelihood):
super(Poisson, self).__init__(gp_link, name='Poisson') super(Poisson, self).__init__(gp_link, name='Poisson')
def _preprocess_values(self,Y): def _conditional_mean(self, f):
return Y """
the expected value of y given a value of f
"""
return self.gp_link.transf(gp)
def pdf_link(self, link_f, y, Y_metadata=None): def pdf_link(self, link_f, y, Y_metadata=None):
""" """