From b856c60d30c7fe33dfd4935b27b4bf42356558d5 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 2 Jul 2013 18:18:11 +0100 Subject: [PATCH] massive changes --- GPy/likelihoods/ep.py | 125 ++++++- GPy/likelihoods/noise_models/__init__.py | 4 + .../noise_models/binomial_likelihood.py | 111 ++++++ .../noise_models/likelihood_functions.py | 349 ++++++++++++++++++ .../noise_models/link_functions.py | 99 +++++ .../noise_models/poisson_likelihood.py | 92 +++++ 6 files changed, 770 insertions(+), 10 deletions(-) create mode 100644 GPy/likelihoods/noise_models/__init__.py create mode 100644 GPy/likelihoods/noise_models/binomial_likelihood.py create mode 100644 GPy/likelihoods/noise_models/likelihood_functions.py create mode 100644 GPy/likelihoods/noise_models/link_functions.py create mode 100644 GPy/likelihoods/noise_models/poisson_likelihood.py diff --git a/GPy/likelihoods/ep.py b/GPy/likelihoods/ep.py index 94f760e9..7e90755e 100644 --- a/GPy/likelihoods/ep.py +++ b/GPy/likelihoods/ep.py @@ -24,9 +24,18 @@ class EP(likelihood): #Initial values - Likelihood approximation parameters: #p(y|f) = t(f|tau_tilde,v_tilde) + #TODO restore self.tau_tilde = np.zeros(self.N) self.v_tilde = np.zeros(self.N) + #_gp = self.LikelihoodFunction.link.transf(self.data) + #_mean = self.LikelihoodFunction._mean(_gp) + #_variance = self.LikelihoodFunction._variance(_gp) + #self.tau_tilde = 1./_variance + #self.tau_tilde[_variance== 0] = 1. + #self.v_tilde = _mean*self.tau_tilde + + #initial values for the GP variables self.Y = np.zeros((self.N,1)) self.covariance_matrix = np.eye(self.N) @@ -38,16 +47,17 @@ class EP(likelihood): self.trYYT = 0. def restart(self): + #FIXME self.tau_tilde = np.zeros(self.N) self.v_tilde = np.zeros(self.N) - 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.V = self.precision * self.Y - self.VVT_factor = self.V - self.trYYT = 0. + #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.V = self.precision * self.Y + #self.VVT_factor = self.V + #self.trYYT = 0. def predictive_values(self,mu,var,full_cov): if full_cov: @@ -78,6 +88,8 @@ class EP(likelihood): self.VVT_factor = self.V self.trYYT = np.trace(self.YYT) + #a = kjkjkjkj + def fit_full(self,K): """ The expectation-propagation algorithm. @@ -117,15 +129,103 @@ class EP(likelihood): self.v_[i] = mu[i]/Sigma[i,i] - self.eta*self.v_tilde[i] #Marginal moments self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.LikelihoodFunction.moments_match(self._transf_data[i],self.tau_[i],self.v_[i]) + + + #DELETE + """ + import pylab as pb + from scipy import stats + import scipy as sp + import link_functions + from constructors import * + + link = link_functions.Log_ex_1() + distribution = poisson(link=link) + gp = np.linspace(-3,50,100) + #distribution = binomial() + #gp = np.linspace(-3,3,100) + + y = self._transf_data[i] + tau_ = self.tau_[i] + v_ = self.v_[i] + sigma2_ = np.sqrt(1./tau_) + mu_ = v_/tau_ + + gaussian = stats.norm.pdf(gp,loc=mu_,scale=np.sqrt(sigma2_)) + non_gaussian = np.array([distribution._mass(gp_i,y) for gp_i in gp]) + prod = np.array([distribution._product(gp_i,y,mu_,np.sqrt(sigma2_)) for gp_i in gp]) + my_Z_hat,my_mu_hat,my_sigma2_hat = distribution.moments_match(y,tau_,v_) + proxy = stats.norm.pdf(gp,loc=my_mu_hat,scale=np.sqrt(my_sigma2_hat)) + + + new_sigma2_tilde = 1./self.tau_tilde[i] + new_mu_tilde = self.v_tilde[i]/self.tau_tilde[i] + new_Z_tilde = self.Z_hat[i]*np.sqrt(2*np.pi)*np.sqrt(sigma2_+new_sigma2_tilde)*np.exp(.5*(mu_-new_mu_tilde)**2/(sigma2_+new_sigma2_tilde)) + bad_gaussian = stats.norm.pdf(gp,self.v_tilde[i]/self.tau_tilde[i],np.sqrt(1./self.tau_tilde[i])) + new_gaussian = stats.norm.pdf(gp,new_mu_tilde,np.sqrt(new_sigma2_tilde))*new_Z_tilde + #new_gaussian = stats.norm.pdf(gp,_mu_tilde,np.sqrt(_sigma2_tilde))*_Z_tilde + + _sigma2_tilde = 1./(1./(my_sigma2_hat) - 1./sigma2_) + _mu_tilde = (my_mu_hat/my_sigma2_hat - mu_/sigma2_)*_sigma2_tilde + _Z_tilde = my_Z_hat*np.sqrt(2*np.pi)*np.sqrt(sigma2_+_sigma2_tilde)*np.exp(.5*(mu_ - _mu_tilde)**2/(sigma2_ + _sigma2_tilde)) + + fig1 = pb.figure(figsize=(15,5)) + ax1 = fig1.add_subplot(131) + ax1.grid(True) + #pb.plot(gp,bad_gaussian,'b--',linewidth=1.5) + #pb.plot(gp,non_gaussian,'b-',linewidth=1.5) + pb.plot(gp,new_gaussian,'r--',linewidth=1.5) + pb.title('Likelihood: $p(y_i|f_i)$',fontsize=22) + + ax2 = fig1.add_subplot(132) + ax2.grid(True) + pb.plot(gp,gaussian,'b-',linewidth=1.5) + pb.title('Cavity distribution: $q_{-i}(f_i)$',fontsize=22) + + ax3 = fig1.add_subplot(133) + ax3.grid(True) + pb.plot(gp,prod,'b--',linewidth=1.5) + + pb.plot(gp,proxy*my_Z_hat,'r-',linewidth=1.5) + + pb.title('Approximation: $\mathcal{N}(f_i|\hat{\mu}_i,\hat{\sigma}_i^2) \hat{Z}_i$',fontsize=22) + pb.legend(('Exact','Approximation'),frameon=False) + + print 'i',i + print 'v/tau _tilde', self.v_tilde[i], self.tau_tilde[i] + print 'v/tau _', self.v_[i], self.tau_[i] + print 'Z/mu/sigma2 _hat', self.Z_hat[i], mu_hat[i], sigma2_hat[i] + pb.plot(gp,new_gaussian*gaussian,'k-') + + a = kj + break + """ + #DELETE + + + + #Site parameters update - Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i]) - Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i]) + Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i]) #FIXME + Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i]) #FIXME self.tau_tilde[i] += Delta_tau self.v_tilde[i] += Delta_v + + #new_tau = self.delta/self.eta*(1./sigma2_hat[i] - self.tau_[i]) + #new_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - self.v_[i]) + #Delta_tau = new_tau - self.tau_tilde[i] + #Delta_v = new_v - self.v_tilde[i] + #self.tau_tilde[i] += Delta_tau + #self.v_tilde[i] += Delta_v + #Posterior distribution parameters update DSYR(Sigma,Sigma[:,i].copy(), -float(Delta_tau/(1.+ Delta_tau*Sigma[i,i]))) mu = np.dot(Sigma,self.v_tilde) self.iterations += 1 + + + + #Sigma recomptutation with Cholesky decompositon Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*K B = np.eye(self.N) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K @@ -138,6 +238,11 @@ class EP(likelihood): self.np1.append(self.tau_tilde.copy()) self.np2.append(self.v_tilde.copy()) + ##DELETE + #pb.vlines(mu[i],0,max(prod)) + #break + #DELETE + return self._compute_GP_variables() def fit_DTC(self, Kmm, Kmn): diff --git a/GPy/likelihoods/noise_models/__init__.py b/GPy/likelihoods/noise_models/__init__.py new file mode 100644 index 00000000..40282402 --- /dev/null +++ b/GPy/likelihoods/noise_models/__init__.py @@ -0,0 +1,4 @@ +import likelihood_functions +import binomial_likelihood +import poisson_likelihood +import link_functions diff --git a/GPy/likelihoods/noise_models/binomial_likelihood.py b/GPy/likelihoods/noise_models/binomial_likelihood.py new file mode 100644 index 00000000..d23dd2f7 --- /dev/null +++ b/GPy/likelihoods/noise_models/binomial_likelihood.py @@ -0,0 +1,111 @@ +# 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 link_functions +from likelihood_functions import NoiseModel + +class Binomial(NoiseModel): + """ + Probit likelihood + Y is expected to take values in {-1,1} + ----- + $$ + L(x) = \\Phi (Y_i*f_i) + $$ + """ + def __init__(self,link=None,analytical_moments=False): + super(Binomial, self).__init__(link,analytical_moments) + + 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.link.d2inv_transf_df2(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 + """ diff --git a/GPy/likelihoods/noise_models/likelihood_functions.py b/GPy/likelihoods/noise_models/likelihood_functions.py new file mode 100644 index 00000000..5376e8e7 --- /dev/null +++ b/GPy/likelihoods/noise_models/likelihood_functions.py @@ -0,0 +1,349 @@ +# 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 GPy.util.plot import gpplot +from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf +import link_functions + + +class NoiseModel(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,analytical_moments=False): + #assert isinstance(link,link_functions.LinkFunction), "link is not a valid LinkFunction."#FIXME + self.link = link + self.analytical_moments = analytical_moments + 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.vstack(pred_mean) + pred_var = np.vstack(pred_var) + q1 = np.vstack(q1) + q3 = np.vstack(q3) + return pred_mean, pred_var, q1, q3 diff --git a/GPy/likelihoods/noise_models/link_functions.py b/GPy/likelihoods/noise_models/link_functions.py new file mode 100644 index 00000000..b0cdcd49 --- /dev/null +++ b/GPy/likelihoods/noise_models/link_functions.py @@ -0,0 +1,99 @@ +# 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 GPy.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(1.+np.exp(f)) + + 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) diff --git a/GPy/likelihoods/noise_models/poisson_likelihood.py b/GPy/likelihoods/noise_models/poisson_likelihood.py new file mode 100644 index 00000000..86a2df2a --- /dev/null +++ b/GPy/likelihoods/noise_models/poisson_likelihood.py @@ -0,0 +1,92 @@ +# 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 GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf +import link_functions +from likelihood_functions import NoiseModel + +class Poisson(NoiseModel): + """ + 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,analytical_moments=False): + #self.discrete = True + #self.support_limits = (0,np.inf) + + #self.analytical_moments = False + super(Poisson, self).__init__(link,analytical_moments) + + def _preprocess_values(self,Y): #TODO + self.scale = .5*Y.max() + self.shift = Y.mean() + return (Y - self.shift)/self.scale + + def _mass(self,gp,obs): + """ + Mass (or density) function + """ + obs = obs*self.scale + self.shift + 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 _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)