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

more files

Browse source
This commit is contained in:
Ricardo 2013-07-01 20:31:16 +01:00
parent 7361d311c1
commit 4054442462
8 changed files with 48 additions and 667 deletions
Download patch file Download diff file Expand all files Collapse all files

6
GPy/examples/regression.py
View file

@ -83,7 +83,7 @@ def coregionalisation_toy2(optim_iters=100):
Y = np.vstack((Y1,Y2)) Y = np.vstack((Y1,Y2))
k1 = GPy.kern.rbf(1) + GPy.kern.bias(1) k1 = GPy.kern.rbf(1) + GPy.kern.bias(1)
k2 = GPy.kern.Coregionalise(2,1) k2 = GPy.kern.coregionalise(2,1)
k = k1.prod(k2,tensor=True) k = k1.prod(k2,tensor=True)
m = GPy.models.GPRegression(X,Y,kernel=k) m = GPy.models.GPRegression(X,Y,kernel=k)
m.constrain_fixed('.*rbf_var',1.) m.constrain_fixed('.*rbf_var',1.)
@ -114,7 +114,7 @@ def coregionalisation_toy(optim_iters=100):
Y = np.vstack((Y1,Y2)) Y = np.vstack((Y1,Y2))
k1 = GPy.kern.rbf(1) k1 = GPy.kern.rbf(1)
k2 = GPy.kern.Coregionalise(2,2) k2 = GPy.kern.coregionalise(2,2)
k = k1.prod(k2,tensor=True) k = k1.prod(k2,tensor=True)
m = GPy.models.GPRegression(X,Y,kernel=k) m = GPy.models.GPRegression(X,Y,kernel=k)
m.constrain_fixed('.*rbf_var',1.) m.constrain_fixed('.*rbf_var',1.)
@ -149,7 +149,7 @@ def coregionalisation_sparse(optim_iters=100):
Z = np.hstack((np.random.rand(num_inducing,1)*8,np.random.randint(0,2,num_inducing)[:,None])) Z = np.hstack((np.random.rand(num_inducing,1)*8,np.random.randint(0,2,num_inducing)[:,None]))
k1 = GPy.kern.rbf(1) k1 = GPy.kern.rbf(1)
k2 = GPy.kern.Coregionalise(2,2) k2 = GPy.kern.coregionalise(2,2)
k = k1.prod(k2,tensor=True) + GPy.kern.white(2,0.001) k = k1.prod(k2,tensor=True) + GPy.kern.white(2,0.001)
m = GPy.models.SparseGPRegression(X,Y,kernel=k,Z=Z) m = GPy.models.SparseGPRegression(X,Y,kernel=k,Z=Z)

5
GPy/likelihoods/__init__.py
View file

@ -1,6 +1,5 @@
from ep import EP from ep import EP
from gaussian import Gaussian from gaussian import Gaussian
# TODO: from Laplace import Laplace # TODO: from Laplace import Laplace
import likelihood_functions as functions from constructors import *
import binomial_likelihood
import poisson_likelihood

123
GPy/likelihoods/binomial_likelihood.py
View file

@ -1,123 +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 pylab as pb
from ..util.plot import gpplot
from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import link_functions
from likelihood_functions import LikelihoodFunction
class Binomial(LikelihoodFunction):
"""
Probit likelihood
Y is expected to take values in {-1,1}
-----
$$
L(x) = \\Phi (Y_i*f_i)
$$
"""
def __init__(self,link=None):
self.discrete = True
self.support_limits = (0,1)
if not link:
link = link_functions.Probit
if isinstance(link,link_functions.Probit):
self.analytical_moments = True
else:
self.analytical_moments = False
super(Binomial, self).__init__(link)
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, 'Binomial 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)
"""
z = data_i*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 + data_i*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)
return Z_hat, mu_hat, sigma2_hat
def _predictive_mean_analytical(self,mu,sigma):
return stats.norm.cdf(mu/np.sqrt(1+sigma**2))
def _mass(self,gp,obs):
#NOTE obs must be in {0,1}
p = self.link.inv_transf(gp)
return p**obs * (1.-p)**(1.-obs)
def _nlog_mass(self,gp,obs):
p = self.link.inv_transf(gp)
return obs*np.log(p) + (1.-obs)*np.log(1-p)
def _dnlog_mass_dgp(self,gp,obs):
p = self.link.inv_transf(gp)
dp = self.link.dinv_transf_df(gp)
return obs/p * dp - (1.-obs)/(1.-p) * dp
def _d2nlog_mass_dgp2(self,gp,obs):
p = self.link.inv_transf(gp)
return (obs/p + (1.-obs)/(1.-p))*self.lind.d2inv_transf_df(gp) + ((1.-obs)/(1.-p)**2-obs/p**2)*self.link.dinv_transf_df(gp)
def _mean(self,gp):
"""
Mass (or density) function
"""
return self.link.inv_transf(gp)
def _dmean_dgp(self,gp):
return self.link.dinv_transf_df(gp)
def _d2mean_dgp2(self,gp):
return self.link.d2inv_transf_df2(gp)
def _variance(self,gp):
"""
Mass (or density) function
"""
p = self.link.inv_transf(gp)
return p*(1-p)
def _dvariance_dgp(self,gp):
return self.link.dinv_transf_df(gp)*(1. - 2.*self.link.inv_transf(gp))
def _d2variance_dgp2(self,gp):
return self.link.d2inv_transf_df2(gp)*(1. - 2.*self.link.inv_transf(gp)) - 2*self.link.dinv_transf_df(gp)**2
"""
def predictive_values(self,mu,var): #TODO remove
mu = mu.flatten()
var = var.flatten()
#mean = stats.norm.cdf(mu/np.sqrt(1+var))
mean = self._predictive_mean_analytical(mu,np.sqrt(var))
norm_025 = [stats.norm.ppf(.025,m,v) for m,v in zip(mu,var)]
norm_975 = [stats.norm.ppf(.975,m,v) for m,v in zip(mu,var)]
#p_025 = stats.norm.cdf(norm_025/np.sqrt(1+var))
#p_975 = stats.norm.cdf(norm_975/np.sqrt(1+var))
p_025 = self._predictive_mean_analytical(norm_025,np.sqrt(var))
p_975 = self._predictive_mean_analytical(norm_975,np.sqrt(var))
return mean[:,None], np.nan*var, p_025[:,None], p_975[:,None] # TODO: var
"""

42
GPy/likelihoods/constructors.py Normal file
View file

@ -0,0 +1,42 @@
# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
from likelihood_functions import LikelihoodFunction
import noise_models
import link_functions
def binomial(link=None):
"""
Construct a binomial likelihood
:param link: a GPy link function
"""
#self.discrete = True
#self.support_limits = (0,1)
if link is None:
link = link_functions.Probit()
else:
assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.'
if isinstance(link,link_functions.Probit):
analytical_moments = True
else:
analytical_moments = False
return noise_models.binomial_likelihood.Binomial(link,analytical_moments)
def poisson(link=None):
"""
Construct a Poisson likelihood
:param link: a GPy link function
"""
if link is None:
link = link_functions.Log_ex_1()
else:
assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.'
#assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.'
analytical_moments = False
return noise_models.poisson_likelihood.Poisson(link,analytical_moments)

348
GPy/likelihoods/likelihood_functions.py
View file

@ -1,348 +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 pylab as pb
from ..util.plot import gpplot
from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import link_functions
class LikelihoodFunction(object):
"""
Likelihood class for doing Expectation propagation
:param Y: observed output (Nx1 numpy.darray)
..Note:: Y values allowed depend on the LikelihoodFunction used
"""
def __init__(self,link):
#assert isinstance(link,link_functions.LinkFunction), "link is not a valid LinkFunction."#FIXME
self.link = link
if self.analytical_moments:
self.moments_match = self._moments_match_analytical
self.predictive_mean = self._predictive_mean_analytical
else:
self.moments_match = self._moments_match_numerical
self.predictive_mean = self._predictive_mean_numerical
def _preprocess_values(self,Y):
"""
In case it is needed, this function assess the output values or makes any pertinent transformation on them.
:param Y: observed output (Nx1 numpy.darray)
"""
return Y
def _product(self,gp,obs,mu,sigma):
"""
Product between the cavity distribution and a likelihood factor.
:param gp: latent variable
:param obs: observed output
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return stats.norm.pdf(gp,loc=mu,scale=sigma) * self._mass(gp,obs)
def _nlog_product_scaled(self,gp,obs,mu,sigma):
"""
Negative log-product between the cavity distribution and a likelihood factor.
..Note:: The constant term in the Gaussian distribution is ignored.
:param gp: latent variable
:param obs: observed output
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return .5*((gp-mu)/sigma)**2 + self._nlog_mass(gp,obs)
def _dnlog_product_dgp(self,gp,obs,mu,sigma):
"""
Derivative wrt latent variable of the log-product between the cavity distribution and a likelihood factor.
:param gp: latent variable
:param obs: observed output
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return (gp - mu)/sigma**2 + self._dnlog_mass_dgp(gp,obs)
def _d2nlog_product_dgp2(self,gp,obs,mu,sigma):
"""
Second derivative wrt latent variable of the log-product between the cavity distribution and a likelihood factor.
:param gp: latent variable
:param obs: observed output
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return 1./sigma**2 + self._d2nlog_mass_dgp2(gp,obs)
def _product_mode(self,obs,mu,sigma):
"""
Newton's CG method to find the mode in _product (cavity x likelihood factor).
:param obs: observed output
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return sp.optimize.fmin_ncg(self._nlog_product_scaled,x0=mu,fprime=self._dnlog_product_dgp,fhess=self._d2nlog_product_dgp2,args=(obs,mu,sigma))
def _moments_match_analytical(self,obs,tau,v):
"""
If available, this function computes the moments analytically.
"""
pass
def _moments_match_numerical(self,obs,tau,v):
"""
Lapace approximation to calculate the moments.
:param obs: observed output
:param tau: cavity distribution 1st natural parameter (precision)
:param v: cavity distribution 2nd natural paramenter (mu*precision)
"""
mu = v/tau
mu_hat = self._product_mode(obs,mu,np.sqrt(1./tau))
sigma2_hat = 1./(tau + self._d2nlog_mass_dgp2(mu_hat,obs))
Z_hat = np.exp(-.5*tau*(mu_hat-mu)**2) * self._mass(mu_hat,obs)*np.sqrt(tau*sigma2_hat)
return Z_hat,mu_hat,sigma2_hat
def _nlog_conditional_mean_scaled(self,gp,mu,sigma):
"""
Negative logarithm of the l.v.'s predictive distribution times the output's mean given the l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
..Note:: This function helps computing E(Y_star) = E(E(Y_star|f_star))
"""
return .5*((gp - mu)/sigma)**2 - np.log(self._mean(gp))
def _dnlog_conditional_mean_dgp(self,gp,mu,sigma):
"""
Derivative of _nlog_conditional_mean_scaled wrt. l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return (gp - mu)/sigma**2 - self._dmean_dgp(gp)/self._mean(gp)
def _d2nlog_conditional_mean_dgp2(self,gp,mu,sigma):
"""
Second derivative of _nlog_conditional_mean_scaled wrt. l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return 1./sigma**2 - self._d2mean_dgp2(gp)/self._mean(gp) + (self._dmean_dgp(gp)/self._mean(gp))**2
def _nlog_exp_conditional_variance_scaled(self,gp,mu,sigma):
"""
Negative logarithm of the l.v.'s predictive distribution times the output's variance given the l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
..Note:: This function helps computing E(V(Y_star|f_star))
"""
return .5*((gp - mu)/sigma)**2 - np.log(self._variance(gp))
def _dnlog_exp_conditional_variance_dgp(self,gp,mu,sigma):
"""
Derivative of _nlog_exp_conditional_variance_scaled wrt. l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return (gp - mu)/sigma**2 - self._dvariance_dgp(gp)/self._variance(gp)
def _d2nlog_exp_conditional_variance_dgp2(self,gp,mu,sigma):
"""
Second derivative of _nlog_exp_conditional_variance_scaled wrt. l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return 1./sigma**2 - self._d2variance_dgp2(gp)/self._variance(gp) + (self._dvariance_dgp(gp)/self._variance(gp))**2
def _nlog_exp_conditional_mean_sq_scaled(self,gp,mu,sigma):
"""
Negative logarithm of the l.v.'s predictive distribution times the output's mean squared given the l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
..Note:: This function helps computing E( E(Y_star|f_star)**2 )
"""
return .5*((gp - mu)/sigma)**2 - 2*np.log(self._mean(gp))
def _dnlog_exp_conditional_mean_sq_dgp(self,gp,mu,sigma):
"""
Derivative of _nlog_exp_conditional_mean_sq_scaled wrt. l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return (gp - mu)/sigma**2 - 2*self._dmean_dgp(gp)/self._mean(gp)
def _d2nlog_exp_conditional_mean_sq_dgp2(self,gp,mu,sigma):
"""
Second derivative of _nlog_exp_conditional_mean_sq_scaled wrt. l.v.
:param gp: latent variable
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
return 1./sigma**2 - 2*( self._d2mean_dgp2(gp)/self._mean(gp) - (self._dmean_dgp(gp)/self._mean(gp))**2 )
def _predictive_mean_analytical(self,mu,sigma):
"""
If available, this function computes the predictive mean analytically.
"""
pass
def _predictive_mean_numerical(self,mu,sigma):
"""
Laplace approximation to the predictive mean: E(Y_star) = E( E(Y_star|f_star) )
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
maximum = sp.optimize.fmin_ncg(self._nlog_conditional_mean_scaled,x0=self._mean(mu),fprime=self._dnlog_conditional_mean_dgp,fhess=self._d2nlog_conditional_mean_dgp2,args=(mu,sigma))
mean = np.exp(-self._nlog_conditional_mean_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_conditional_mean_dgp2(maximum,mu,sigma))*sigma)
"""
pb.figure()
x = np.array([mu + step*sigma for step in np.linspace(-7,7,100)])
f = np.array([np.exp(-self._nlog_conditional_mean_scaled(xi,mu,sigma))/np.sqrt(2*np.pi*sigma**2) for xi in x])
pb.plot(x,f,'b-')
sigma2 = 1./self._d2nlog_conditional_mean_dgp2(maximum,mu,sigma)
f2 = np.exp(-.5*(x-maximum)**2/sigma2)/np.sqrt(2*np.pi*sigma2)
k = np.exp(-self._nlog_conditional_mean_scaled(maximum,mu,sigma))*np.sqrt(sigma2)/np.sqrt(sigma**2)
pb.plot(x,f2*mean,'r-')
pb.vlines(maximum,0,f.max())
"""
return mean
def _predictive_mean_sq(self,mu,sigma):
"""
Laplace approximation to the predictive mean squared: E(Y_star**2) = E( E(Y_star|f_star)**2 )
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
"""
maximum = sp.optimize.fmin_ncg(self._nlog_exp_conditional_mean_sq_scaled,x0=self._mean(mu),fprime=self._dnlog_exp_conditional_mean_sq_dgp,fhess=self._d2nlog_exp_conditional_mean_sq_dgp2,args=(mu,sigma))
mean_squared = np.exp(-self._nlog_exp_conditional_mean_sq_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_exp_conditional_mean_sq_dgp2(maximum,mu,sigma))*sigma)
return mean_squared
def predictive_variance(self,mu,sigma,predictive_mean=None):
"""
Laplace approximation to the predictive variance: V(Y_star) = E( V(Y_star|f_star) ) + V( E(Y_star|f_star) )
:param mu: cavity distribution mean
:param sigma: cavity distribution standard deviation
:predictive_mean: output's predictive mean, if None _predictive_mean function will be called.
"""
# E( V(Y_star|f_star) )
maximum = sp.optimize.fmin_ncg(self._nlog_exp_conditional_variance_scaled,x0=self._variance(mu),fprime=self._dnlog_exp_conditional_variance_dgp,fhess=self._d2nlog_exp_conditional_variance_dgp2,args=(mu,sigma))
exp_var = np.exp(-self._nlog_exp_conditional_variance_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_exp_conditional_variance_dgp2(maximum,mu,sigma))*sigma)
"""
pb.figure()
x = np.array([mu + step*sigma for step in np.linspace(-7,7,100)])
f = np.array([np.exp(-self._nlog_exp_conditional_variance_scaled(xi,mu,sigma))/np.sqrt(2*np.pi*sigma**2) for xi in x])
pb.plot(x,f,'b-')
sigma2 = 1./self._d2nlog_exp_conditional_variance_dgp2(maximum,mu,sigma)
f2 = np.exp(-.5*(x-maximum)**2/sigma2)/np.sqrt(2*np.pi*sigma2)
k = np.exp(-self._nlog_exp_conditional_variance_scaled(maximum,mu,sigma))*np.sqrt(sigma2)/np.sqrt(sigma**2)
pb.plot(x,f2*exp_var,'r--')
pb.vlines(maximum,0,f.max())
"""
#V( E(Y_star|f_star) ) = E( E(Y_star|f_star)**2 ) - E( E(Y_star|f_star)**2 )
exp_exp2 = self._predictive_mean_sq(mu,sigma)
if predictive_mean is None:
predictive_mean = self.predictive_mean(mu,sigma)
var_exp = exp_exp2 - predictive_mean**2
return exp_var + var_exp
def _nlog_joint_predictive_scaled(self,x,mu,sigma):
"""
Negative logarithm of the joint predictive distribution (latent variable and output).
:param x: tuple (latent variable,output)
:param mu: latent variable's predictive mean
:param sigma: latent variable's predictive standard deviation
"""
return self._nlog_product_scaled(x[0],x[1],mu,sigma)
def _gradient_nlog_joint_predictive(self,x,mu,sigma):
"""
Gradient of _nlog_joint_predictive_scaled.
:param x: tuple (latent variable,output)
:param mu: latent variable's predictive mean
:param sigma: latent variable's predictive standard deviation
..Note: Only avilable when the output is continuous
"""
assert not self.discrete, "Gradient not available for discrete outputs."
return np.array((self._dnlog_product_dgp(gp=x[0],obs=x[1],mu=mu,sigma=sigma),self._dnlog_mass_dobs(obs=x[1],gp=x[0])))
def _hessian_nlog_joint_predictive(self,x,mu,sigma):
"""
Hessian of _nlog_joint_predictive_scaled.
:param x: tuple (latent variable,output)
:param mu: latent variable's predictive mean
:param sigma: latent variable's predictive standard deviation
..Note: Only avilable when the output is continuous
"""
assert not self.discrete, "Hessian not available for discrete outputs."
cross_derivative = self._d2nlog_mass_dcross(gp=x[0],obs=x[1])
return np.array((self._d2nlog_product_dgp2(gp=x[0],obs=x[1],mu=mu,sigma=sigma),cross_derivative,cross_derivative,self._d2nlog_mass_dobs2(obs=x[1],gp=x[0]))).reshape(2,2)
def _joint_predictive_mode(self,mu,sigma):
"""
Negative logarithm of the joint predictive distribution (latent variable and output).
:param x: tuple (latent variable,output)
:param mu: latent variable's predictive mean
:param sigma: latent variable's predictive standard deviation
"""
return sp.optimize.fmin_ncg(self._nlog_joint_predictive_scaled,x0=(mu,self.link.inv_transf(mu)),fprime=self._gradient_nlog_joint_predictive,fhess=self._hessian_nlog_joint_predictive,args=(mu,sigma))
def predictive_values(self,mu,var,sample=True,sample_size=5000):
"""
Compute mean, variance and conficence interval (percentiles 5 and 95) of the prediction
:param mu: mean of the latent variable
:param var: variance of the latent variable
"""
if isinstance(mu,float) or isinstance(mu,int):
mu = [mu]
var = [var]
pred_mean = []
pred_var = []
q1 = []
q3 = []
for m,s in zip(mu,np.sqrt(var)):
pred_mean.append(self.predictive_mean(m,s))
pred_var.append(self.predictive_variance(m,s,pred_mean[-1]))
q1.append(self.predictive_mean(stats.norm.ppf(.025,m,s**2),s))
q3.append(self.predictive_mean(stats.norm.ppf(.975,m,s**2),s))
pred_mean = np.array(pred_mean)[:,None]
pred_var = np.array(pred_var)[:,None]
q1 = np.array(q1)[:,None]
q3 = np.array(q3)[:,None]
return pred_mean, pred_var, q1, q3

100
GPy/likelihoods/link_functions.py
View file

@ -1,100 +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
import scipy as sp
import pylab as pb
from ..util.plot import gpplot
from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf,inv_std_norm_cdf
class LinkFunction(object):
"""
Link function class for doing non-Gaussian likelihoods approximation
:param Y: observed output (Nx1 numpy.darray)
..Note:: Y values allowed depend on the likelihood_function used
"""
def __init__(self):
pass
class Identity(LinkFunction):
"""
$$
g(f) = f
$$
"""
def transf(self,mu):
return mu
def inv_transf(self,f):
return f
def dinv_transf_df(self,f):
return 1.
def d2inv_transf_df2(self,f):
return 0
class Probit(LinkFunction):
"""
$$
g(f) = \\Phi^{-1} (mu)
$$
"""
def transf(self,mu):
return inv_std_norm_cdf(mu)
def inv_transf(self,f):
return std_norm_cdf(f)
def dinv_transf_df(self,f):
return std_norm_pdf(f)
def d2inv_transf_df2(self,f):
return -f * std_norm_pdf(f)
class Log(LinkFunction):
"""
$$
g(f) = \log(\mu)
$$
"""
def transf(self,mu):
return np.log(mu)
def inv_transf(self,f):
return np.exp(f)
def dinv_transf_df(self,f):
return np.exp(f)
def d2inv_transf_df2(self,f):
return np.exp(f)
class Log_ex_1(LinkFunction):
"""
$$
g(f) = \log(\exp(\mu) - 1)
$$
"""
def transf(self,mu):
"""
function: output space -> latent space
"""
return np.log(np.exp(mu) - 1)
def inv_transf(self,f):
"""
function: latent space -> output space
"""
return np.log(np.exp(f)+1)
def dinv_transf_df(self,f):
return np.exp(f)/(1.+np.exp(f))
def d2inv_transf_df2(self,f):
aux = np.exp(f)/(1.+np.exp(f))
return aux*(1.-aux)

89
GPy/likelihoods/poisson_likelihood.py
View file

@ -1,89 +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 pylab as pb
from ..util.plot import gpplot
from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import link_functions
from likelihood_functions import LikelihoodFunction
class Poisson(LikelihoodFunction):
"""
Poisson likelihood
Y is expected to take values in {0,1,2,...}
-----
$$
L(x) = \exp(\lambda) * \lambda**Y_i / Y_i!
$$
"""
def __init__(self,link=None):
self.discrete = True
self.support_limits = (0,np.inf)
self.analytical_moments = False
super(Poisson, self).__init__(link)
def _mass(self,gp,obs):
"""
Mass (or density) function
"""
return stats.poisson.pmf(obs,self.link.inv_transf(gp))
def _nlog_mass(self,gp,obs):
"""
Negative logarithm of the un-normalized distribution: factors that are not a function of gp are omitted
"""
return self.link.inv_transf(gp) - obs * np.log(self.link.inv_transf(gp)) + np.log(special.gamma(obs+1))
#def _preprocess_values(self,Y): #TODO
def _dnlog_mass_dgp(self,gp,obs):
return self.link.dinv_transf_df(gp) * (1. - obs/self.link.inv_transf(gp))
def _d2nlog_mass_dgp2(self,gp,obs):
d2_df = self.link.d2inv_transf_df2(gp)
inv_transf = self.link.inv_transf(gp)
return obs * ((self.link.dinv_transf_df(gp)/inv_transf)**2 - d2_df/inv_transf) + d2_df
def _dnlog_mass_dobs(self,obs,gp): #TODO not needed
return special.psi(obs+1) - np.log(self.link.inv_transf(gp))
def _d2nlog_mass_dobs2(self,obs,gp=None): #TODO not needed
return special.polygamma(1,obs)
def _d2nlog_mass_dcross(self,obs,gp): #TODO not needed
return -self.link.dinv_transf_df(gp)/self.link.inv_transf(gp)
def _mean(self,gp):
"""
Mass (or density) function
"""
return self.link.inv_transf(gp)
#def _variance(self,gp):
# return self.link.inv_transf(gp)
def _dmean_dgp(self,gp):
return self.link.dinv_transf_df(gp)
def _d2mean_dgp2(self,gp):
return self.link.d2inv_transf_df2(gp)
def _variance(self,gp):
"""
Mass (or density) function
"""
return self.link.inv_transf(gp)
#def _variance(self,gp):
# return self.link.inv_transf(gp)
def _dvariance_dgp(self,gp):
return self.link.dinv_transf_df(gp)
def _d2variance_dgp2(self,gp):
return self.link.d2inv_transf_df2(gp)

2
GPy/models/gp_classification.py
View file

@ -32,7 +32,7 @@ class GPClassification(GP):
if likelihood is None: if likelihood is None:
#distribution = GPy.likelihoods.binomial_likelihood.Binomial(link=link) #distribution = GPy.likelihoods.binomial_likelihood.Binomial(link=link)
distribution = likelihoods.binomial_likelihood.Binomial() distribution = likelihoods.binomial()
likelihood = likelihoods.EP(Y, distribution) likelihood = likelihoods.EP(Y, distribution)
elif Y is not None: elif Y is not None:
if not all(Y.flatten() == likelihood.data.flatten()): if not all(Y.flatten() == likelihood.data.flatten()):