apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-04-28 06:16:24 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Bug fixed in numerical approx. to the predictive variance.

Browse source
This commit is contained in:
Ricardo 2013-11-11 08:39:58 +00:00
parent c6014b0f17
commit 604e60d5cf
2 changed files with 11 additions and 12 deletions
Download patch file Download diff file Expand all files Collapse all files

2
GPy/likelihoods/noise_models/gp_transformations.py
View file

@ -78,9 +78,11 @@ class Probit(GPTransformation):
return std_norm_pdf(f)
def d2transf_df2(self,f):
#FIXME
return -f * std_norm_pdf(f)
def d3transf_df3(self,f):
#FIXME
f2 = f**2
return -(1/(np.sqrt(2*np.pi)))*np.exp(-0.5*(f2))*(1-f2)

21
GPy/likelihoods/noise_models/noise_distributions.py
View file

@ -99,31 +99,29 @@ class NoiseDistribution(object):
:param tau: cavity distribution 1st natural parameter (precision)
:param v: cavity distribution 2nd natural paramenter (mu*precision)
"""
#Compute first integral for zeroth moment
#Compute first integral for zeroth moment.
#NOTE constant np.sqrt(2*pi/tau) added at the end of the function
mu = v/tau
def int_1(f):
return self.pdf(f, obs)*np.exp(-0.5*tau*np.square(mu-f))
z, accuracy = quad(int_1, -np.inf, np.inf)
#z /= np.sqrt(2*np.pi/tau)
z_scaled, accuracy = quad(int_1, -np.inf, np.inf)
#Compute second integral for first moment
def int_2(f):
return f*self.pdf(f, obs)*np.exp(-0.5*tau*np.square(mu-f))
mean, accuracy = quad(int_2, -np.inf, np.inf)
#mean /= np.sqrt(2*np.pi/tau)
mean /= z
mean /= z_scaled
#Compute integral for variance
def int_3(f):
return (f**2)*self.pdf(f, obs)*np.exp(-0.5*tau*np.square(mu-f))
Ef2, accuracy = quad(int_3, -np.inf, np.inf)
#Ef2 /= np.sqrt(2*np.pi/tau)
Ef2 /= z
Ef2 /= z_scaled
variance = Ef2 - mean**2
#Add constant to the zeroth moment
#NOTE: this constant is not needed in the other moments because it cancells out.
z /= np.sqrt(2*np.pi/tau)
z = z_scaled/np.sqrt(2*np.pi/tau)
return z, mean, variance
@ -185,18 +183,17 @@ class NoiseDistribution(object):
#V( E(Y_star|f_star) ) = E( E(Y_star|f_star)**2 ) - E( E(Y_star|f_star) )**2
#E( E(Y_star|f_star)**2 )
#E( E(Y_star|f_star) )**2
if predictive_mean is None:
predictive_mean = self.predictive_mean(mu,variance)
predictive_mean_sq = predictive_mean**2
#E( E(Y_star|f_star)**2 )
def int_pred_mean_sq(f,m,v,predictive_mean_sq):
return predictive_mean_sq*np.exp(-(0.5/v)*np.square(f - m))
return self._mean(f)**2*np.exp(-(0.5/v)*np.square(f - m))
scaled_exp_exp2 = [quad(int_pred_mean_sq, -np.inf, np.inf,args=(mj,s2j,pm2j))[0] for mj,s2j,pm2j in zip(mu,variance,predictive_mean_sq)]
exp_exp2 = np.array(scaled_exp_exp2)[:,None] / normalizer
#E( E(Y_star|f_star) )**2
var_exp = exp_exp2 - predictive_mean_sq
# V(Y_star) = E( V(Y_star|f_star) ) + V( E(Y_star|f_star) )