mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-07 02:52:40 +02:00
massive restructuting to make the EP likelihoods work consistently
This commit is contained in:
James Hensman
parent
ea0802d938
commit
a6851cf63d
2 changed files with 67 additions and 86 deletions
|
|
@ -1,313 +0,0 @@
|
|||
import numpy as np
|
||||
import random
|
||||
import pylab as pb #TODO erase me
|
||||
from scipy import stats, linalg
|
||||
from .likelihoods import likelihood
|
||||
from ..core import model
|
||||
from ..util.linalg import pdinv,mdot,jitchol
|
||||
from ..util.plot import gpplot
|
||||
from .. import kern
|
||||
|
||||
class EP:
|
||||
def __init__(self,covariance,likelihood,Kmn=None,Knn_diag=None,epsilon=1e-3,power_ep=[1.,1.]):
|
||||
"""
|
||||
Expectation Propagation
|
||||
|
||||
Arguments
|
||||
---------
|
||||
X : input observations
|
||||
likelihood : Output's likelihood (likelihood class)
|
||||
kernel : a GPy kernel (kern class)
|
||||
inducing : Either an array specifying the inducing points location or a sacalar defining their number. None value for using a non-sparse model is used.
|
||||
power_ep : Power-EP parameters (eta,delta) - 2x1 numpy array (floats)
|
||||
epsilon : Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float)
|
||||
"""
|
||||
self.likelihood = likelihood
|
||||
assert covariance.shape[0] == covariance.shape[1]
|
||||
if Kmn is not None:
|
||||
self.Kmm = covariance
|
||||
self.Kmn = Kmn
|
||||
self.M = self.Kmn.shape[0]
|
||||
self.N = self.Kmn.shape[1]
|
||||
assert self.M < self.N, 'The number of inducing inputs must be smaller than the number of observations'
|
||||
else:
|
||||
self.K = covariance
|
||||
self.N = self.K.shape[0]
|
||||
if Knn_diag is not None:
|
||||
self.Knn_diag = Knn_diag
|
||||
assert len(Knn_diag) == self.N, 'Knn_diagonal has size different from N'
|
||||
|
||||
self.epsilon = epsilon
|
||||
self.eta, self.delta = power_ep
|
||||
self.jitter = 1e-12
|
||||
|
||||
"""
|
||||
Initial values - Likelihood approximation parameters:
|
||||
p(y|f) = t(f|tau_tilde,v_tilde)
|
||||
"""
|
||||
self.tau_tilde = np.zeros(self.N)
|
||||
self.v_tilde = np.zeros(self.N)
|
||||
|
||||
def _compute_GP_variables(self):
|
||||
#Variables to be called from GP
|
||||
mu_tilde = self.v_tilde/self.tau_tilde #When calling EP, this variable is used instead of Y in the GP model
|
||||
sigma_sum = 1./self.tau_ + 1./self.tau_tilde
|
||||
mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2
|
||||
Z_ep = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant
|
||||
return self.tau_tilde[:,None], mu_tilde[:,None], Z_ep
|
||||
|
||||
class Full(EP):
|
||||
def fit_EP(self):
|
||||
"""
|
||||
The expectation-propagation algorithm.
|
||||
For nomenclature see Rasmussen & Williams 2006.
|
||||
"""
|
||||
#Prior distribution parameters: p(f|X) = N(f|0,K)
|
||||
#self.K = self.kernel.K(self.X,self.X)
|
||||
|
||||
#Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma)
|
||||
self.mu=np.zeros(self.N)
|
||||
self.Sigma=self.K.copy()
|
||||
|
||||
"""
|
||||
Initial values - Cavity distribution parameters:
|
||||
q_(f|mu_,sigma2_) = Product{q_i(f|mu_i,sigma2_i)}
|
||||
sigma_ = 1./tau_
|
||||
mu_ = v_/tau_
|
||||
"""
|
||||
self.tau_ = np.empty(self.N,dtype=float)
|
||||
self.v_ = np.empty(self.N,dtype=float)
|
||||
|
||||
#Initial values - Marginal moments
|
||||
z = np.empty(self.N,dtype=float)
|
||||
self.Z_hat = np.empty(self.N,dtype=float)
|
||||
phi = np.empty(self.N,dtype=float)
|
||||
mu_hat = np.empty(self.N,dtype=float)
|
||||
sigma2_hat = np.empty(self.N,dtype=float)
|
||||
|
||||
#Approximation
|
||||
epsilon_np1 = self.epsilon + 1.
|
||||
epsilon_np2 = self.epsilon + 1.
|
||||
self.iterations = 0
|
||||
self.np1 = [self.tau_tilde.copy()]
|
||||
self.np2 = [self.v_tilde.copy()]
|
||||
while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon:
|
||||
update_order = np.arange(self.N)
|
||||
random.shuffle(update_order)
|
||||
for i in update_order:
|
||||
#Cavity distribution parameters
|
||||
self.tau_[i] = 1./self.Sigma[i,i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = self.mu[i]/self.Sigma[i,i] - self.eta*self.v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.likelihood.moments_match(i,self.tau_[i],self.v_[i])
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./self.Sigma[i,i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - self.mu[i]/self.Sigma[i,i])
|
||||
self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
|
||||
self.v_tilde[i] = self.v_tilde[i] + Delta_v
|
||||
#Posterior distribution parameters update
|
||||
si=self.Sigma[:,i].reshape(self.N,1)
|
||||
self.Sigma = self.Sigma - Delta_tau/(1.+ Delta_tau*self.Sigma[i,i])*np.dot(si,si.T)
|
||||
self.mu = np.dot(self.Sigma,self.v_tilde)
|
||||
self.iterations += 1
|
||||
#Sigma recomptutation with Cholesky decompositon
|
||||
Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*(self.K)
|
||||
B = np.eye(self.N) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K
|
||||
L = jitchol(B)
|
||||
V,info = linalg.flapack.dtrtrs(L,Sroot_tilde_K,lower=1)
|
||||
self.Sigma = self.K - np.dot(V.T,V)
|
||||
self.mu = np.dot(self.Sigma,self.v_tilde)
|
||||
epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.N
|
||||
epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N
|
||||
self.np1.append(self.tau_tilde.copy())
|
||||
self.np2.append(self.v_tilde.copy())
|
||||
|
||||
return self._compute_GP_variables()
|
||||
|
||||
class DTC(EP):
|
||||
def fit_EP(self):
|
||||
"""
|
||||
The expectation-propagation algorithm with sparse pseudo-input.
|
||||
For nomenclature see ... 2013.
|
||||
"""
|
||||
|
||||
"""
|
||||
Prior approximation parameters:
|
||||
q(f|X) = int_{df}{N(f|KfuKuu_invu,diag(Kff-Qff)*N(u|0,Kuu)} = N(f|0,Sigma0)
|
||||
Sigma0 = Qnn = Knm*Kmmi*Kmn
|
||||
"""
|
||||
self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm)
|
||||
self.KmnKnm = np.dot(self.Kmn, self.Kmn.T)
|
||||
self.KmmiKmn = np.dot(self.Kmmi,self.Kmn)
|
||||
self.Qnn_diag = np.sum(self.Kmn*self.KmmiKmn,-2)
|
||||
self.LLT0 = self.Kmm.copy()
|
||||
|
||||
"""
|
||||
Posterior approximation: q(f|y) = N(f| mu, Sigma)
|
||||
Sigma = Diag + P*R.T*R*P.T + K
|
||||
mu = w + P*gamma
|
||||
"""
|
||||
self.mu = np.zeros(self.N)
|
||||
self.LLT = self.Kmm.copy()
|
||||
self.Sigma_diag = self.Qnn_diag.copy()
|
||||
|
||||
"""
|
||||
Initial values - Cavity distribution parameters:
|
||||
q_(g|mu_,sigma2_) = Product{q_i(g|mu_i,sigma2_i)}
|
||||
sigma_ = 1./tau_
|
||||
mu_ = v_/tau_
|
||||
"""
|
||||
self.tau_ = np.empty(self.N,dtype=float)
|
||||
self.v_ = np.empty(self.N,dtype=float)
|
||||
|
||||
#Initial values - Marginal moments
|
||||
z = np.empty(self.N,dtype=float)
|
||||
self.Z_hat = np.empty(self.N,dtype=float)
|
||||
phi = np.empty(self.N,dtype=float)
|
||||
mu_hat = np.empty(self.N,dtype=float)
|
||||
sigma2_hat = np.empty(self.N,dtype=float)
|
||||
|
||||
#Approximation
|
||||
epsilon_np1 = 1
|
||||
epsilon_np2 = 1
|
||||
self.iterations = 0
|
||||
self.np1 = [self.tau_tilde.copy()]
|
||||
self.np2 = [self.v_tilde.copy()]
|
||||
while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon:
|
||||
update_order = np.arange(self.N)
|
||||
random.shuffle(update_order)
|
||||
for i in update_order:
|
||||
#Cavity distribution parameters
|
||||
self.tau_[i] = 1./self.Sigma_diag[i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = self.mu[i]/self.Sigma_diag[i] - self.eta*self.v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.likelihood.moments_match(i,self.tau_[i],self.v_[i])
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./self.Sigma_diag[i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - self.mu[i]/self.Sigma_diag[i])
|
||||
self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
|
||||
self.v_tilde[i] = self.v_tilde[i] + Delta_v
|
||||
#Posterior distribution parameters update
|
||||
self.LLT = self.LLT + np.outer(self.Kmn[:,i],self.Kmn[:,i])*Delta_tau
|
||||
L = jitchol(self.LLT)
|
||||
V,info = linalg.flapack.dtrtrs(L,self.Kmn,lower=1)
|
||||
self.Sigma_diag = np.sum(V*V,-2)
|
||||
si = np.sum(V.T*V[:,i],-1)
|
||||
self.mu = self.mu + (Delta_v-Delta_tau*self.mu[i])*si
|
||||
self.iterations += 1
|
||||
#Sigma recomputation with Cholesky decompositon
|
||||
self.LLT0 = self.LLT0 + np.dot(self.Kmn*self.tau_tilde[None,:],self.Kmn.T)
|
||||
self.L = jitchol(self.LLT)
|
||||
V,info = linalg.flapack.dtrtrs(L,self.Kmn,lower=1)
|
||||
V2,info = linalg.flapack.dtrtrs(L.T,V,lower=0)
|
||||
self.Sigma_diag = np.sum(V*V,-2)
|
||||
Knmv_tilde = np.dot(self.Kmn,self.v_tilde)
|
||||
self.mu = np.dot(V2.T,Knmv_tilde)
|
||||
epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.N
|
||||
epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N
|
||||
self.np1.append(self.tau_tilde.copy())
|
||||
self.np2.append(self.v_tilde.copy())
|
||||
|
||||
return self._compute_GP_variables()
|
||||
|
||||
class FITC(EP):
|
||||
def fit_EP(self):
|
||||
"""
|
||||
The expectation-propagation algorithm with sparse pseudo-input.
|
||||
For nomenclature see Naish-Guzman and Holden, 2008.
|
||||
"""
|
||||
|
||||
"""
|
||||
Prior approximation parameters:
|
||||
q(f|X) = int_{df}{N(f|KfuKuu_invu,diag(Kff-Qff)*N(u|0,Kuu)} = N(f|0,Sigma0)
|
||||
Sigma0 = diag(Knn-Qnn) + Qnn, Qnn = Knm*Kmmi*Kmn
|
||||
"""
|
||||
self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm)
|
||||
self.P0 = self.Kmn.T
|
||||
self.KmnKnm = np.dot(self.P0.T, self.P0)
|
||||
self.KmmiKmn = np.dot(self.Kmmi,self.P0.T)
|
||||
self.Qnn_diag = np.sum(self.P0.T*self.KmmiKmn,-2)
|
||||
self.Diag0 = self.Knn_diag - self.Qnn_diag
|
||||
self.R0 = jitchol(self.Kmmi).T
|
||||
|
||||
"""
|
||||
Posterior approximation: q(f|y) = N(f| mu, Sigma)
|
||||
Sigma = Diag + P*R.T*R*P.T + K
|
||||
mu = w + P*gamma
|
||||
"""
|
||||
self.w = np.zeros(self.N)
|
||||
self.gamma = np.zeros(self.M)
|
||||
self.mu = np.zeros(self.N)
|
||||
self.P = self.P0.copy()
|
||||
self.R = self.R0.copy()
|
||||
self.Diag = self.Diag0.copy()
|
||||
self.Sigma_diag = self.Knn_diag
|
||||
|
||||
"""
|
||||
Initial values - Cavity distribution parameters:
|
||||
q_(g|mu_,sigma2_) = Product{q_i(g|mu_i,sigma2_i)}
|
||||
sigma_ = 1./tau_
|
||||
mu_ = v_/tau_
|
||||
"""
|
||||
self.tau_ = np.empty(self.N,dtype=float)
|
||||
self.v_ = np.empty(self.N,dtype=float)
|
||||
|
||||
#Initial values - Marginal moments
|
||||
z = np.empty(self.N,dtype=float)
|
||||
self.Z_hat = np.empty(self.N,dtype=float)
|
||||
phi = np.empty(self.N,dtype=float)
|
||||
mu_hat = np.empty(self.N,dtype=float)
|
||||
sigma2_hat = np.empty(self.N,dtype=float)
|
||||
|
||||
#Approximation
|
||||
epsilon_np1 = 1
|
||||
epsilon_np2 = 1
|
||||
self.iterations = 0
|
||||
self.np1 = [self.tau_tilde.copy()]
|
||||
self.np2 = [self.v_tilde.copy()]
|
||||
while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon:
|
||||
update_order = np.arange(self.N)
|
||||
random.shuffle(update_order)
|
||||
for i in update_order:
|
||||
#Cavity distribution parameters
|
||||
self.tau_[i] = 1./self.Sigma_diag[i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = self.mu[i]/self.Sigma_diag[i] - self.eta*self.v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.likelihood.moments_match(i,self.tau_[i],self.v_[i])
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./self.Sigma_diag[i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - self.mu[i]/self.Sigma_diag[i])
|
||||
self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
|
||||
self.v_tilde[i] = self.v_tilde[i] + Delta_v
|
||||
#Posterior distribution parameters update
|
||||
dtd1 = Delta_tau*self.Diag[i] + 1.
|
||||
dii = self.Diag[i]
|
||||
self.Diag[i] = dii - (Delta_tau * dii**2.)/dtd1
|
||||
pi_ = self.P[i,:].reshape(1,self.M)
|
||||
self.P[i,:] = pi_ - (Delta_tau*dii)/dtd1 * pi_
|
||||
Rp_i = np.dot(self.R,pi_.T)
|
||||
RTR = np.dot(self.R.T,np.dot(np.eye(self.M) - Delta_tau/(1.+Delta_tau*self.Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),self.R))
|
||||
self.R = jitchol(RTR).T
|
||||
self.w[i] = self.w[i] + (Delta_v - Delta_tau*self.w[i])*dii/dtd1
|
||||
self.gamma = self.gamma + (Delta_v - Delta_tau*self.mu[i])*np.dot(RTR,self.P[i,:].T)
|
||||
self.RPT = np.dot(self.R,self.P.T)
|
||||
self.Sigma_diag = self.Diag + np.sum(self.RPT.T*self.RPT.T,-1)
|
||||
self.mu = self.w + np.dot(self.P,self.gamma)
|
||||
self.iterations += 1
|
||||
#Sigma recomptutation with Cholesky decompositon
|
||||
self.Diag = self.Diag0/(1.+ self.Diag0 * self.tau_tilde)
|
||||
self.P = (self.Diag / self.Diag0)[:,None] * self.P0
|
||||
self.RPT0 = np.dot(self.R0,self.P0.T)
|
||||
L = jitchol(np.eye(self.M) + np.dot(self.RPT0,(1./self.Diag0 - self.Diag/(self.Diag0**2))[:,None]*self.RPT0.T))
|
||||
self.R,info = linalg.flapack.dtrtrs(L,self.R0,lower=1)
|
||||
self.RPT = np.dot(self.R,self.P.T)
|
||||
self.Sigma_diag = self.Diag + np.sum(self.RPT.T*self.RPT.T,-1)
|
||||
self.w = self.Diag * self.v_tilde
|
||||
self.gamma = np.dot(self.R.T, np.dot(self.RPT,self.v_tilde))
|
||||
self.mu = self.w + np.dot(self.P,self.gamma)
|
||||
epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.N
|
||||
epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N
|
||||
self.np1.append(self.tau_tilde.copy())
|
||||
self.np2.append(self.v_tilde.copy())
|
||||
|
||||
return self._compute_GP_variables()
|
||||
|
|
@ -1,229 +0,0 @@
|
|||
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||
# 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
|
||||
|
||||
class likelihood:
|
||||
"""
|
||||
Likelihood class for doing Expectation propagation
|
||||
|
||||
:param Y: observed output (Nx1 numpy.darray)
|
||||
..Note:: Y values allowed depend on the likelihood used
|
||||
"""
|
||||
def __init__(self,Y,location=0,scale=1):
|
||||
self.Y = Y
|
||||
self.N = self.Y.shape[0]
|
||||
self.location = location
|
||||
self.scale = scale
|
||||
|
||||
def plot2D(self,X,X_new,F_new,U=None):
|
||||
"""
|
||||
Predictive distribution of the fitted GP model for 2-dimensional inputs
|
||||
|
||||
:param X_new: The points at which to make a prediction
|
||||
:param Mean_new: mean values at X_new
|
||||
:param Var_new: variance values at X_new
|
||||
:param X_u: input points used to train the model
|
||||
:param Mean_u: mean values at X_u
|
||||
:param Var_new: variance values at X_u
|
||||
"""
|
||||
N,D = X_new.shape
|
||||
assert D == 2, 'Number of dimensions must be 2'
|
||||
n = np.sqrt(N)
|
||||
x1min = X_new[:,0].min()
|
||||
x1max = X_new[:,0].max()
|
||||
x2min = X_new[:,1].min()
|
||||
x2max = X_new[:,1].max()
|
||||
pb.imshow(F_new.reshape(n,n),extent=(x1min,x1max,x2max,x2min),vmin=0,vmax=1)
|
||||
pb.colorbar()
|
||||
C1 = np.arange(self.N)[self.Y.flatten()==1]
|
||||
C2 = np.arange(self.N)[self.Y.flatten()==-1]
|
||||
[pb.plot(X[i,0],X[i,1],'ro') for i in C1]
|
||||
[pb.plot(X[i,0],X[i,1],'bo') for i in C2]
|
||||
pb.xlim(x1min,x1max)
|
||||
pb.ylim(x2min,x2max)
|
||||
if U is not None:
|
||||
[pb.plot(a,b,'wo') for a,b in U]
|
||||
|
||||
class probit(likelihood):
|
||||
"""
|
||||
Probit likelihood
|
||||
Y is expected to take values in {-1,1}
|
||||
-----
|
||||
$$
|
||||
L(x) = \\Phi (Y_i*f_i)
|
||||
$$
|
||||
"""
|
||||
def __init__(self,Y,location=0,scale=1):
|
||||
assert np.sum(np.abs(Y)-1) == 0, "Output values must be either -1 or 1"
|
||||
likelihood.__init__(self,Y,location,scale)
|
||||
|
||||
def moments_match(self,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 = self.Y[i]*v_i/np.sqrt(tau_i**2 + tau_i)
|
||||
Z_hat = stats.norm.cdf(z)
|
||||
phi = stats.norm.pdf(z)
|
||||
mu_hat = v_i/tau_i + self.Y[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(self,mu,var):
|
||||
mu = mu.flatten()
|
||||
var = var.flatten()
|
||||
return stats.norm.cdf(mu/np.sqrt(1+var))
|
||||
|
||||
def predictive_var(self,mu,var):
|
||||
p=self.predictive_mean(mu,var)
|
||||
return p*(1-p)
|
||||
|
||||
def _log_likelihood_gradients():
|
||||
raise NotImplementedError
|
||||
|
||||
def plot(self,X,mu,var,phi,X_obs,Z=None,samples=0):
|
||||
assert X_obs.shape[1] == 1, 'Number of dimensions must be 1'
|
||||
phi_var = self.predictive_var(mu,var)
|
||||
gpplot(X,phi,phi_var)
|
||||
if samples:
|
||||
phi_samples = np.vstack([np.random.binomial(1,phi.flatten()) for s in range(samples)])
|
||||
pb.plot(X,phi_samples.T,'x', alpha = 0.4, c='#3465a4' )
|
||||
pb.plot(X_obs,(self.Y+1)/2,'kx',mew=1.5)
|
||||
if Z is not None:
|
||||
pb.plot(Z,Z*0+.5,'r|',mew=1.5,markersize=12)
|
||||
pb.ylim(-0.2,1.2)
|
||||
|
||||
class poisson(likelihood):
|
||||
"""
|
||||
Poisson likelihood
|
||||
Y is expected to take values in {0,1,2,...}
|
||||
-----
|
||||
$$
|
||||
L(x) = \exp(\lambda) * \lambda**Y_i / Y_i!
|
||||
$$
|
||||
"""
|
||||
def __init__(self,Y,location=0,scale=1):
|
||||
assert len(Y[Y<0]) == 0, "Output cannot have negative values"
|
||||
likelihood.__init__(self,Y,location,scale)
|
||||
|
||||
def moments_match(self,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)
|
||||
"""
|
||||
mu = v_i/tau_i
|
||||
sigma = np.sqrt(1./tau_i)
|
||||
def poisson_norm(f):
|
||||
"""
|
||||
Product of the likelihood and the cavity distribution
|
||||
"""
|
||||
pdf_norm_f = stats.norm.pdf(f,loc=mu,scale=sigma)
|
||||
rate = np.exp( (f*self.scale)+self.location)
|
||||
poisson = stats.poisson.pmf(float(self.Y[i]),rate)
|
||||
return pdf_norm_f*poisson
|
||||
|
||||
def log_pnm(f):
|
||||
"""
|
||||
Log of poisson_norm
|
||||
"""
|
||||
return -(-.5*(f-mu)**2/sigma**2 - np.exp( (f*self.scale)+self.location) + ( (f*self.scale)+self.location)*self.Y[i])
|
||||
|
||||
"""
|
||||
Golden Search and Simpson's Rule
|
||||
--------------------------------
|
||||
Simpson's Rule is used to calculate the moments mumerically, it needs a grid of points as input.
|
||||
Golden Search is used to find the mode in the poisson_norm distribution and define around it the grid for Simpson's Rule
|
||||
"""
|
||||
#TODO golden search & simpson's rule can be defined in the general likelihood class, rather than in each specific case.
|
||||
|
||||
#Golden search
|
||||
golden_A = -1 if self.Y[i] == 0 else np.array([np.log(self.Y[i]),mu]).min() #Lower limit
|
||||
golden_B = np.array([np.log(self.Y[i]),mu]).max() #Upper limit
|
||||
golden_A = (golden_A - self.location)/self.scale
|
||||
golden_B = (golden_B - self.location)/self.scale
|
||||
opt = sp.optimize.golden(log_pnm,brack=(golden_A,golden_B)) #Better to work with log_pnm than with poisson_norm
|
||||
|
||||
# Simpson's approximation
|
||||
width = 3./np.log(max(self.Y[i],2))
|
||||
A = opt - width #Lower limit
|
||||
B = opt + width #Upper limit
|
||||
K = 10*int(np.log(max(self.Y[i],150))) #Number of points in the grid, we DON'T want K to be the same number for every case
|
||||
h = (B-A)/K # length of the intervals
|
||||
grid_x = np.hstack([np.linspace(opt-width,opt,K/2+1)[1:-1], np.linspace(opt,opt+width,K/2+1)]) # grid of points (X axis)
|
||||
x = np.hstack([A,B,grid_x[range(1,K,2)],grid_x[range(2,K-1,2)]]) # grid_x rearranged, just to make Simpson's algorithm easier
|
||||
zeroth = np.hstack([poisson_norm(A),poisson_norm(B),[4*poisson_norm(f) for f in grid_x[range(1,K,2)]],[2*poisson_norm(f) for f in grid_x[range(2,K-1,2)]]]) # grid of points (Y axis) rearranged like x
|
||||
first = zeroth*x
|
||||
second = first*x
|
||||
Z_hat = sum(zeroth)*h/3 # Zero-th moment
|
||||
mu_hat = sum(first)*h/(3*Z_hat) # First moment
|
||||
m2 = sum(second)*h/(3*Z_hat) # Second moment
|
||||
sigma2_hat = m2 - mu_hat**2 # Second central moment
|
||||
return float(Z_hat), float(mu_hat), float(sigma2_hat)
|
||||
|
||||
def predictive_mean(self,mu,var):
|
||||
return np.exp(mu*self.scale + self.location)
|
||||
|
||||
def predictive_var(self,mu,var):
|
||||
return predictive_mean(mu,var)
|
||||
|
||||
def _log_likelihood_gradients():
|
||||
raise NotImplementedError
|
||||
|
||||
def plot(self,X,mu,var,phi,X_obs,Z=None,samples=0):
|
||||
assert X_obs.shape[1] == 1, 'Number of dimensions must be 1'
|
||||
gpplot(X,phi,phi.flatten())
|
||||
pb.plot(X_obs,self.Y,'kx',mew=1.5)
|
||||
if samples:
|
||||
phi_samples = np.vstack([np.random.poisson(phi.flatten(),phi.size) for s in range(samples)])
|
||||
pb.plot(X,phi_samples.T, alpha = 0.4, c='#3465a4', linewidth = 0.8)
|
||||
if Z is not None:
|
||||
pb.plot(Z,Z*0+pb.ylim()[0],'k|',mew=1.5,markersize=12)
|
||||
|
||||
class gaussian(likelihood):
|
||||
"""
|
||||
Gaussian likelihood
|
||||
Y is expected to take values in (-inf,inf)
|
||||
"""
|
||||
self.variance = variance
|
||||
self._data = Y
|
||||
self.
|
||||
def moments_match(self,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)
|
||||
"""
|
||||
mu = v_i/tau_i
|
||||
sigma = np.sqrt(1./tau_i)
|
||||
s = 1. if self.Y[i] == 0 else 1./self.Y[i]
|
||||
sigma2_hat = 1./(1./sigma**2 + 1./s**2)
|
||||
mu_hat = sigma2_hat*(mu/sigma**2 + self.Y[i]/s**2)
|
||||
Z_hat = 1./np.sqrt(2*np.pi) * 1./np.sqrt(sigma**2+s**2) * np.exp(-.5*(mu-self.Y[i])**2/(sigma**2 + s**2))
|
||||
return Z_hat, mu_hat, sigma2_hat
|
||||
|
||||
def plot1Db(self,X,X_new,F_new,U=None):
|
||||
assert X.shape[1] == 1, 'Number of dimensions must be 1'
|
||||
gpplot(X_new,F_new,np.zeros(X_new.shape[0]))
|
||||
pb.plot(X,self.Y,'kx',mew=1.5)
|
||||
if U is not None:
|
||||
pb.plot(U,np.ones(U.shape[0])*self.Y.min()*.8,'r|',mew=1.5,markersize=12)
|
||||
|
||||
def _log_likelihood_gradients():
|
||||
raise NotImplementedError
|
||||
else:
|
||||
var = var[:,None] * np.square(self._Ystd)
|
||||
|
||||
Loading…
Add table
Add a link
Reference in a new issue