Merge pull request #534 from adhaka/EP_refactor

Ep refactor
This commit is contained in:
Zhenwen Dai 2017-09-04 10:39:30 +01:00 committed by GitHub
commit 44a232ff6a
19 changed files with 517 additions and 76 deletions

View file

@ -16,8 +16,9 @@ addons:
env:
- PYTHON_VERSION=2.7
#- PYTHON_VERSION=3.3
- PYTHON_VERSION=3.4
#- PYTHON_VERSION=3.4
- PYTHON_VERSION=3.5
- PYTHON_VERSION=3.6
before_install:
- wget https://github.com/mzwiessele/travis_scripts/raw/master/download_miniconda.sh

View file

@ -1,5 +1,97 @@
# Changelog
## v1.7.6 (2017-06-19)
### Fix
* Appveyor not uploading to testpypi for now. [mzwiessele]
### Other
* Bump version: 1.7.5 → 1.7.6. [mzwiessele]
## v1.7.5 (2017-06-19)
### Fix
* Splitting forecast tests into 3 to circumvent 10 minute stop of travis. [mzwiessele]
### Other
* Bump version: 1.7.4 → 1.7.5. [mzwiessele]
## v1.7.4 (2017-06-19)
### Fix
* Paramz version for parallel optimization fix. [mzwiessele]
### Other
* Bump version: 1.7.3 → 1.7.4. [mzwiessele]
## v1.7.3 (2017-06-19)
### Fix
* Appveyor build failing. [mzwiessele]
### Other
* Bump version: 1.7.2 → 1.7.3. [mzwiessele]
## v1.7.2 (2017-06-17)
### Fix
* Appveyor build python 3.6. [mzwiessele]
### Other
* Bump version: 1.7.1 → 1.7.2. [mzwiessele]
## v1.7.1 (2017-06-17)
### Fix
* Appveyor build python 3.6. [mzwiessele]
### Other
* Bump version: 1.7.0 → 1.7.1. [mzwiessele]
## v1.7.0 (2017-06-17)
### Fix
* Support for 3.5 and higher now that 3.6 is out. [mzwiessele]
### Other
* Bump version: 1.6.3 → 1.7.0. [mzwiessele]
## v1.6.3 (2017-06-17)
### Other
* Bump version: 1.6.2 → 1.6.3. [mzwiessele]
* Merge pull request #504 from rmcantin/devel. [Max Zwiessele]
* Fix python 2-3 compatibility. [Ruben Martinez-Cantin]
* Merge pull request #511 from dirmeier/devel. [Max Zwiessele]
* Added LICENSE file to MANIFEST.in. [dirmeier]
## v1.6.2 (2017-04-12)
### Fix

View file

@ -1 +1 @@
__version__ = "1.6.2"
__version__ = "1.7.7"

View file

@ -6,6 +6,7 @@ from paramz import ObsAr
from . import ExactGaussianInference, VarDTC
from ...util import diag
from .posterior import PosteriorEP as Posterior
from ...likelihoods import Gaussian
log_2_pi = np.log(2*np.pi)
@ -174,18 +175,18 @@ class EP(EPBase, ExactGaussianInference):
if self.ep_mode=="nested":
#Force EP at each step of the optimization
self._ep_approximation = None
post_params, ga_approx, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
post_params, ga_approx, cav_params, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
elif self.ep_mode=="alternated":
if getattr(self, '_ep_approximation', None) is None:
#if we don't yet have the results of runnign EP, run EP and store the computed factors in self._ep_approximation
post_params, ga_approx, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
post_params, ga_approx, cav_params, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
else:
#if we've already run EP, just use the existing approximation stored in self._ep_approximation
post_params, ga_approx, log_Z_tilde = self._ep_approximation
post_params, ga_approx, cav_params, log_Z_tilde = self._ep_approximation
else:
raise ValueError("ep_mode value not valid")
return self._inference(K, ga_approx, likelihood, Y_metadata=Y_metadata, Z_tilde=log_Z_tilde)
return self._inference(Y, K, ga_approx, cav_params, likelihood, Y_metadata=Y_metadata, Z_tilde=log_Z_tilde)
def expectation_propagation(self, K, Y, likelihood, Y_metadata):
@ -220,7 +221,7 @@ class EP(EPBase, ExactGaussianInference):
# This terms cancel with the coreresponding terms in the marginal loglikelihood
log_Z_tilde = self._log_Z_tilde(marg_moments, ga_approx, cav_params)
# - 0.5*np.log(tau_tilde) + 0.5*(v_tilde*v_tilde*1./tau_tilde)
return (post_params, ga_approx, log_Z_tilde)
return (post_params, ga_approx, cav_params, log_Z_tilde)
def _init_approximations(self, K, num_data):
#initial values - Gaussian factors
@ -280,7 +281,7 @@ class EP(EPBase, ExactGaussianInference):
return log_marginal, post_params
def _inference(self, K, ga_approx, likelihood, Z_tilde, Y_metadata=None):
def _inference(self, Y, K, ga_approx, cav_params, likelihood, Z_tilde, Y_metadata=None):
log_marginal, post_params = self._ep_marginal(K, ga_approx, Z_tilde)
tau_tilde_root = np.sqrt(ga_approx.tau)
@ -293,8 +294,7 @@ class EP(EPBase, ExactGaussianInference):
symmetrify(Wi) #(K + Sigma^(\tilde))^(-1)
dL_dK = 0.5 * (tdot(alpha) - Wi)
dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK), Y_metadata)
dL_dthetaL = likelihood.ep_gradients(Y, cav_params.tau, cav_params.v, np.diag(dL_dK), Y_metadata=Y_metadata, quad_mode='gh')
return Posterior(woodbury_inv=Wi, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL, 'dL_dm':alpha}

View file

@ -66,7 +66,14 @@ class Binomial(Likelihood):
np.testing.assert_array_equal(N.shape, y.shape)
nchoosey = special.gammaln(N+1) - special.gammaln(y+1) - special.gammaln(N-y+1)
return nchoosey + y*np.log(inv_link_f) + (N-y)*np.log(1.-inv_link_f)
Ny = N-y
t1 = np.zeros(y.shape)
t2 = np.zeros(y.shape)
t1[y>0] = y[y>0]*np.log(inv_link_f[y>0])
t2[Ny>0] = Ny[Ny>0]*np.log(1.-inv_link_f[Ny>0])
return nchoosey + t1 + t2
def dlogpdf_dlink(self, inv_link_f, y, Y_metadata=None):
"""
@ -86,7 +93,13 @@ class Binomial(Likelihood):
N = Y_metadata['trials']
np.testing.assert_array_equal(N.shape, y.shape)
return y/inv_link_f - (N-y)/(1.-inv_link_f)
Ny = N-y
t1 = np.zeros(y.shape)
t2 = np.zeros(y.shape)
t1[y>0] = y[y>0]/inv_link_f[y>0]
t2[Ny>0] = (Ny[Ny>0])/(1.-inv_link_f[Ny>0])
return t1 - t2
def d2logpdf_dlink2(self, inv_link_f, y, Y_metadata=None):
"""
@ -111,7 +124,13 @@ class Binomial(Likelihood):
"""
N = Y_metadata['trials']
np.testing.assert_array_equal(N.shape, y.shape)
return -y/np.square(inv_link_f) - (N-y)/np.square(1.-inv_link_f)
Ny = N-y
t1 = np.zeros(y.shape)
t2 = np.zeros(y.shape)
t1[y>0] = -y[y>0]/np.square(inv_link_f[y>0])
t2[Ny>0] = -(Ny[Ny>0])/np.square(1.-inv_link_f[Ny>0])
return t1+t2
def d3logpdf_dlink3(self, inv_link_f, y, Y_metadata=None):
"""
@ -135,8 +154,14 @@ class Binomial(Likelihood):
N = Y_metadata['trials']
np.testing.assert_array_equal(N.shape, y.shape)
inv_link_f2 = np.square(inv_link_f)
return 2*y/inv_link_f**3 - 2*(N-y)/(1.-inv_link_f)**3
#inv_link_f2 = np.square(inv_link_f) #TODO Remove. Why is this here?
Ny = N-y
t1 = np.zeros(y.shape)
t2 = np.zeros(y.shape)
t1[y>0] = 2*y[y>0]/inv_link_f[y>0]**3
t2[Ny>0] = - 2*(Ny[Ny>0])/(1.-inv_link_f[Ny>0])**3
return t1 + t2
def samples(self, gp, Y_metadata=None, **kw):
"""

View file

@ -57,7 +57,10 @@ class Gaussian(Likelihood):
def update_gradients(self, grad):
self.variance.gradient = grad
def exact_inference_gradients(self, dL_dKdiag,Y_metadata=None):
def ep_gradients(self, Y, cav_tau, cav_v, dL_dKdiag, Y_metadata=None, quad_mode='gk', boost_grad=1.):
return self.exact_inference_gradients(dL_dKdiag)
def exact_inference_gradients(self, dL_dKdiag, Y_metadata=None):
return dL_dKdiag.sum()
def _preprocess_values(self, Y):

View file

@ -6,8 +6,12 @@ from scipy import stats,special
import scipy as sp
from . import link_functions
from ..util.misc import chain_1, chain_2, chain_3, blockify_dhess_dtheta, blockify_third, blockify_hessian, safe_exp
from ..util.quad_integrate import quadgk_int
from scipy.integrate import quad
from functools import partial
import warnings
from ..core.parameterization import Parameterized
class Likelihood(Parameterized):
@ -223,6 +227,91 @@ class Likelihood(Parameterized):
self.__gh_points = np.polynomial.hermite.hermgauss(T)
return self.__gh_points
def ep_gradients(self, Y, cav_tau, cav_v, dL_dKdiag, Y_metadata=None, quad_mode='gk', boost_grad=1.):
if self.size > 0:
shape = Y.shape
tau,v,Y = cav_tau.flatten(), cav_v.flatten(),Y.flatten()
mu = v/tau
sigma2 = 1./tau
# assert Y.shape == v.shape
dlik_dtheta = np.empty((self.size, Y.shape[0]))
# for j in range(self.size):
Y_metadata_list = []
for index in range(len(Y)):
Y_metadata_i = {}
if Y_metadata is not None:
for key in Y_metadata.keys():
Y_metadata_i[key] = Y_metadata[key][index,:]
Y_metadata_list.append(Y_metadata_i)
if quad_mode == 'gk':
f = partial(self.integrate_gk)
quads = zip(*map(f, Y.flatten(), mu.flatten(), np.sqrt(sigma2.flatten()), Y_metadata_list))
quads = np.vstack(quads)
quads.reshape(self.size, shape[0], shape[1])
elif quad_mode == 'gh':
f = partial(self.integrate_gh)
quads = zip(*map(f, Y.flatten(), mu.flatten(), np.sqrt(sigma2.flatten())))
quads = np.hstack(quads)
quads = quads.T
else:
raise Exception("no other quadrature mode available")
# do a gaussian-hermite integration
dL_dtheta_avg = boost_grad * np.nanmean(quads, axis=1)
dL_dtheta = boost_grad * np.nansum(quads, axis=1)
# dL_dtheta = boost_grad * np.nansum(dlik_dtheta, axis=1)
else:
dL_dtheta = np.zeros(self.num_params)
return dL_dtheta
def integrate_gk(self, Y, mu, sigma, Y_metadata_i=None):
# gaussian-kronrod integration.
fmin = -np.inf
fmax = np.inf
SQRT_2PI = np.sqrt(2.*np.pi)
def generate_integral(f):
a = np.exp(self.logpdf_link(f, Y, Y_metadata_i)) * np.exp(-0.5 * np.square((f - mu) / sigma)) / (
SQRT_2PI * sigma)
fn1 = a * self.dlogpdf_dtheta(f, Y, Y_metadata_i)
fn = fn1
return fn
dF_dtheta_i = quadgk_int(generate_integral, fmin=fmin, fmax=fmax)
return dF_dtheta_i
def integrate_gh(self, Y, mu, sigma, Y_metadata_i=None, gh_points=None):
# gaussian-hermite quadrature.
# "calculate site derivatives E_f{d logp(y_i|f_i)/da} where a is a likelihood parameter
# and the expectation is over the exact marginal posterior, which is not gaussian- and is
# unnormalised product of the cavity distribution(a Gaussian) and the exact likelihood term.
#
# calculate the expectation wrt the approximate marginal posterior, which should be approximately the same.
# . This term is needed for evaluating the
# gradients of the marginal likelihood estimate Z_EP wrt likelihood parameters."
# "writing it explicitly "
# use them for gaussian-hermite quadrature
SQRT_2PI = np.sqrt(2.*np.pi)
if gh_points is None:
gh_x, gh_w = self._gh_points(32)
else:
gh_x, gh_w = gh_points
X = gh_x[None,:]*np.sqrt(2.)*sigma + mu
# Here X is a grid vector of possible fi values, while Y is just a single value which will be broadcasted.
a = np.exp(self.logpdf_link(X, Y, Y_metadata_i))
a = a.repeat(self.num_params,0)
b = self.dlogpdf_dtheta(X, Y, Y_metadata_i)
old_shape = b.shape
fn = np.array([i*j for i,j in zip(a.flatten(), b.flatten())])
fn = fn.reshape(old_shape)
dF_dtheta_i = np.dot(fn, gh_w)/np.sqrt(np.pi)
return dF_dtheta_i
def variational_expectations(self, Y, m, v, gh_points=None, Y_metadata=None):
"""
Use Gauss-Hermite Quadrature to compute

View file

@ -28,10 +28,10 @@ class TestObservationModels(unittest.TestCase):
self.Y_noisy = self.Y.copy()
self.Y_verynoisy = self.Y.copy()
self.Y_noisy[75:80] += 1.3
self.Y_noisy[75] += 1.3
self.init_var = 0.3
self.deg_free = 5.
self.init_var = 0.15
self.deg_free = 4.
censored = np.zeros_like(self.Y)
random_inds = np.random.choice(self.N, int(self.N / 2), replace=True)
censored[random_inds] = 1
@ -83,7 +83,7 @@ class TestObservationModels(unittest.TestCase):
# taking laplace predictions as the ground truth
probs_mean_lap, probs_var_lap = m1.predict(self.X)
probs_mean_ep_alt, probs_var_ep_alt = m2.predict(self.X)
probs_mean_ep_nested, probs_var_ep_nested = m2.predict(self.X)
probs_mean_ep_nested, probs_var_ep_nested = m3.predict(self.X)
# for simple single dimension data , marginal likelihood for laplace and EP approximations should not be so far apart.
self.assertAlmostEqual(m1.log_likelihood(), m2.log_likelihood(),delta=1)
@ -107,12 +107,12 @@ class TestObservationModels(unittest.TestCase):
ep_inf_nested = GPy.inference.latent_function_inference.EP(ep_mode='nested')
ep_inf_frac = GPy.inference.latent_function_inference.EP(ep_mode='nested', eta=0.7)
m1 = GPy.core.GP(self.X, self.Y_noisy.copy(), kernel=self.kernel1, likelihood=studentT.copy(), inference_method=laplace_inf)
m1 = GPy.core.GP(self.X.copy(), self.Y_noisy.copy(), kernel=self.kernel1.copy(), likelihood=studentT.copy(), inference_method=laplace_inf)
# optimize
m1['.*white'].constrain_fixed(1e-5)
m1.randomize()
m2 = GPy.core.GP(self.X, self.Y_noisy.copy(), kernel=self.kernel1, likelihood=studentT.copy(), inference_method=ep_inf_alt)
m2 = GPy.core.GP(self.X.copy(), self.Y_noisy.copy(), kernel=self.kernel1.copy(), likelihood=studentT.copy(), inference_method=ep_inf_alt)
m2['.*white'].constrain_fixed(1e-5)
# m2.constrain_bounded('.*t_scale2', 0.001, 10)
m2.randomize()
@ -124,20 +124,22 @@ class TestObservationModels(unittest.TestCase):
optimizer='bfgs'
m1.optimize(optimizer=optimizer,max_iters=400)
m2.optimize(optimizer=optimizer, max_iters=500)
m2.optimize(optimizer=optimizer, max_iters=400)
# m3.optimize(optimizer=optimizer, max_iters=500)
self.assertAlmostEqual(m1.log_likelihood(), m2.log_likelihood(),delta=200)
self.assertAlmostEqual(m1.log_likelihood(), m2.log_likelihood(),delta=10)
# self.assertAlmostEqual(m1.log_likelihood(), m3.log_likelihood(), 3)
preds_mean_lap, preds_var_lap = m1.predict(self.X)
preds_mean_alt, preds_var_alt = m2.predict(self.X)
# preds_mean_nested, preds_var_nested = m3.predict(self.X)
rmse_lap = self.rmse(preds_mean_lap, self.Y_noisy)
rmse_alt = self.rmse(preds_mean_alt, self.Y_noisy)
rmse_lap = self.rmse(preds_mean_lap, self.Y)
rmse_alt = self.rmse(preds_mean_alt, self.Y)
# rmse_nested = self.rmse(preds_mean_nested, self.Y_noisy)
if rmse_alt > rmse_alt:
self.assertAlmostEqual(rmse_lap, rmse_alt, delta=1.)
if rmse_alt > rmse_lap:
self.assertAlmostEqual(rmse_lap, rmse_alt, delta=1.5)
# m3.optimize(optimizer=optimizer, max_iters=500)

View file

@ -306,11 +306,7 @@ class StateSpaceKernelsTests(np.testing.TestCase):
gp_kernel=gp_kernel,
mean_compare_decimal=2, var_compare_decimal=2)
def test_forecast(self,):
"""
Test time-series forecasting.
"""
def test_forecast_regular(self,):
# Generate data ->
np.random.seed(339) # seed the random number generator
#import pdb; pdb.set_trace()
@ -334,37 +330,102 @@ class StateSpaceKernelsTests(np.testing.TestCase):
#import pdb; pdb.set_trace()
def get_new_kernels():
periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
return ss_kernel, gp_kernel
ss_kernel, gp_kernel = get_new_kernels()
self.run_for_model(X_train, Y_train, ss_kernel, kalman_filter_type = 'regular',
use_cython=False, optimize_max_iters=30, check_gradients=True,
predict_X=X_test,
gp_kernel=gp_kernel,
mean_compare_decimal=2, var_compare_decimal=2)
def test_forecast_svd(self,):
# Generate data ->
np.random.seed(339) # seed the random number generator
#import pdb; pdb.set_trace()
(X,Y) = generate_sine_data(x_points=None, sin_period=5.0, sin_ampl=5.0, noise_var=2.0,
plot = False, points_num=100, x_interval = (0, 40), random=True)
(X1,Y1) = generate_linear_data(x_points=X, tangent=1.0, add_term=20.0, noise_var=0.0,
plot = False, points_num=100, x_interval = (0, 40), random=True)
Y = Y + Y1
X_train = X[X <= 20]
Y_train = Y[X <= 20]
X_test = X[X > 20]
Y_test = Y[X > 20]
X.shape = (X.shape[0],1); Y.shape = (Y.shape[0],1)
X_train.shape = (X_train.shape[0],1); Y_train.shape = (Y_train.shape[0],1)
X_test.shape = (X_test.shape[0],1); Y_test.shape = (Y_test.shape[0],1)
# Generate data <-
#import pdb; pdb.set_trace()
periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
ss_kernel, gp_kernel = get_new_kernels()
self.run_for_model(X_train, Y_train, ss_kernel, kalman_filter_type = 'svd',
use_cython=False, optimize_max_iters=30, check_gradients=False,
predict_X=X_test,
gp_kernel=gp_kernel,
mean_compare_decimal=2, var_compare_decimal=2)
ss_kernel, gp_kernel = get_new_kernels()
def test_forecast_svd_cython(self,):
# Generate data ->
np.random.seed(339) # seed the random number generator
#import pdb; pdb.set_trace()
(X,Y) = generate_sine_data(x_points=None, sin_period=5.0, sin_ampl=5.0, noise_var=2.0,
plot = False, points_num=100, x_interval = (0, 40), random=True)
(X1,Y1) = generate_linear_data(x_points=X, tangent=1.0, add_term=20.0, noise_var=0.0,
plot = False, points_num=100, x_interval = (0, 40), random=True)
Y = Y + Y1
X_train = X[X <= 20]
Y_train = Y[X <= 20]
X_test = X[X > 20]
Y_test = Y[X > 20]
X.shape = (X.shape[0],1); Y.shape = (Y.shape[0],1)
X_train.shape = (X_train.shape[0],1); Y_train.shape = (Y_train.shape[0],1)
X_test.shape = (X_test.shape[0],1); Y_test.shape = (Y_test.shape[0],1)
# Generate data <-
#import pdb; pdb.set_trace()
periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
self.run_for_model(X_train, Y_train, ss_kernel, kalman_filter_type = 'svd',
use_cython=True, optimize_max_iters=30, check_gradients=False,
predict_X=X_test,

View file

@ -64,13 +64,13 @@ class InferenceGPEP(unittest.TestCase):
def genNoisyData(self):
np.random.seed(1)
X = np.random.rand(100,1)
self.real_std = 0.2
self.real_std = 0.1
noise = np.random.randn(*X[:, 0].shape)*self.real_std
Y = (np.sin(X[:, 0]*2*np.pi) + noise)[:, None]
self.f = np.random.rand(X.shape[0],1)
Y_extra_noisy = Y.copy()
Y_extra_noisy[50:53] += 4.
Y_extra_noisy[80:83] -= 2.
Y_extra_noisy[50] += 4.
# Y_extra_noisy[80:83] -= 2.
return X, Y, Y_extra_noisy
def test_inference_EP(self):
@ -85,10 +85,11 @@ class InferenceGPEP(unittest.TestCase):
inference_method=inf,
likelihood=lik)
K = self.model.kern.K(X)
post_params, ga_approx, log_Z_tilde = self.model.inference_method.expectation_propagation(K, ObsAr(Y), lik, None)
post_params, ga_approx, cav_params, log_Z_tilde = self.model.inference_method.expectation_propagation(K, ObsAr(Y), lik, None)
mu_tilde = ga_approx.v / ga_approx.tau.astype(float)
p, m, d = self.model.inference_method._inference(K, ga_approx, lik, Y_metadata=None, Z_tilde=log_Z_tilde)
p, m, d = self.model.inference_method._inference(Y, K, ga_approx, cav_params, lik, Y_metadata=None, Z_tilde=log_Z_tilde)
p0, m0, d0 = super(GPy.inference.latent_function_inference.expectation_propagation.EP, inf).inference(k, X,lik ,mu_tilde[:,None], mean_function=None, variance=1./ga_approx.tau, K=K, Z_tilde=log_Z_tilde + np.sum(- 0.5*np.log(ga_approx.tau) + 0.5*(ga_approx.v*ga_approx.v*1./ga_approx.tau)))
assert (np.sum(np.array([m - m0,
@ -109,19 +110,19 @@ class InferenceGPEP(unittest.TestCase):
def test_inference_EP_non_classification(self):
from paramz import ObsAr
X, Y, Y_extra_noisy = self.genNoisyData()
deg_freedom = 5
init_noise_var = 0.4
deg_freedom = 5.
init_noise_var = 0.08
lik_studentT = GPy.likelihoods.StudentT(deg_free=deg_freedom, sigma2=init_noise_var)
# like_gaussian_noise = GPy.likelihoods.MixedNoise()
k = GPy.kern.RBF(1, variance=2., lengthscale=1.1)
ep_inf_alt = GPy.inference.latent_function_inference.expectation_propagation.EP(max_iters=100, delta=0.5)
ep_inf_nested = GPy.inference.latent_function_inference.expectation_propagation.EP(ep_mode='nested', max_iters=100, delta=0.5)
ep_inf_alt = GPy.inference.latent_function_inference.expectation_propagation.EP(max_iters=4, delta=0.5)
# ep_inf_nested = GPy.inference.latent_function_inference.expectation_propagation.EP(ep_mode='nested', max_iters=100, delta=0.5)
m = GPy.core.GP(X=X,Y=Y_extra_noisy,kernel=k,likelihood=lik_studentT,inference_method=ep_inf_alt)
K = m.kern.K(X)
post_params, ga_approx, log_Z_tilde = m.inference_method.expectation_propagation(K, ObsAr(Y_extra_noisy), lik_studentT, None)
post_params, ga_approx, cav_params, log_Z_tilde = m.inference_method.expectation_propagation(K, ObsAr(Y_extra_noisy), lik_studentT, None)
mu_tilde = ga_approx.v / ga_approx.tau.astype(float)
p, m, d = m.inference_method._inference(K, ga_approx, lik_studentT, Y_metadata=None, Z_tilde=log_Z_tilde)
p, m, d = m.inference_method._inference(Y_extra_noisy, K, ga_approx, cav_params, lik_studentT, Y_metadata=None, Z_tilde=log_Z_tilde)
p0, m0, d0 = super(GPy.inference.latent_function_inference.expectation_propagation.EP, ep_inf_alt).inference(k, X,lik_studentT ,mu_tilde[:,None], mean_function=None, variance=1./ga_approx.tau, K=K, Z_tilde=log_Z_tilde + np.sum(- 0.5*np.log(ga_approx.tau) + 0.5*(ga_approx.v*ga_approx.v*1./ga_approx.tau)))
assert (np.sum(np.array([m - m0,

View file

@ -0,0 +1,39 @@
from __future__ import print_function, division
import numpy as np
import GPy
import warnings
from ..util.quad_integrate import quadgk_int, quadvgk
class QuadTests(np.testing.TestCase):
"""
test file for checking implementation of gaussian-kronrod quadrature.
we will take a function which can be integrated analytically and check if quadgk result is similar or not!
through this file we can test how numerically accurate quadrature implementation in native numpy or manual code is.
"""
def setUp(self):
pass
def test_infinite_quad(self):
def f(x):
return np.exp(-0.5*x**2)*np.power(x,np.arange(3)[:,None])
quad_int_val = quadgk_int(f)
real_val = np.sqrt(np.pi * 2)
np.testing.assert_almost_equal(real_val, quad_int_val[0], decimal=7)
def test_finite_quad(self):
def f2(x):
return x**2
quad_int_val = quadvgk(f2, 1.,2.)
real_val = 7/3.
np.testing.assert_almost_equal(real_val, quad_int_val, decimal=5)
if __name__ == '__main__':
def f(x):
return np.exp(-0.5 * x ** 2) * np.power(x, np.arange(3)[:, None])
quad_int_val = quadgk_int(f)
real_val = np.sqrt(np.pi*2)
np.testing.assert_almost_equal(real_val, quad_int_val[0], decimal=7)
print(quadgk_int(f))

View file

@ -17,3 +17,4 @@ from . import multioutput
from . import parallel
from . import functions
from . import cluster_with_offset
from . import quad_integrate

View file

@ -206,7 +206,10 @@ def authorize_download(dataset_name=None):
def download_data(dataset_name=None):
"""Check with the user that the are happy with terms and conditions for the data set, then download it."""
import itertools
try:
from itertools import zip_longest
except ImportError:
from itertools import izip_longest as zip_longest
dr = data_resources[dataset_name]
if not authorize_download(dataset_name):
@ -220,8 +223,8 @@ def download_data(dataset_name=None):
if 'suffices' in dr: zip_urls += (dr['suffices'], )
else: zip_urls += ([],)
for url, files, save_names, suffices in itertools.zip_longest(*zip_urls, fillvalue=[]):
for f, save_name, suffix in itertools.zip_longest(files, save_names, suffices, fillvalue=None):
for url, files, save_names, suffices in zip_longest(*zip_urls, fillvalue=[]):
for f, save_name, suffix in zip_longest(files, save_names, suffices, fillvalue=None):
download_url(os.path.join(url,f), dataset_name, save_name, suffix=suffix)
return True

119
GPy/util/quad_integrate.py Normal file
View file

@ -0,0 +1,119 @@
"""
The file for utilities related to integration by quadrature methods
- will contain implementation for gaussian-kronrod integration.
"""
import numpy as np
def getSubs(Subs, XK, NK=1):
M = (Subs[1, :] - Subs[0, :]) / 2
C = (Subs[1, :] + Subs[0, :]) / 2
I = XK[:, None] * M + np.ones((NK, 1)) * C
# A = [Subs(1,:); I]
A = np.vstack((Subs[0, :], I))
# B = [I;Subs(2,:)]
B = np.vstack((I, Subs[1, :]))
# Subs = [reshape(A, 1, []);
A = A.flatten()
# reshape(B, 1, [])];
B = B.flatten()
Subs = np.vstack((A,B))
# Subs = np.concatenate((A, B), axis=0)
return Subs
def quadvgk(feval, fmin, fmax, tol1=1e-5, tol2=1e-5):
"""
numpy implementation makes use of the code here: http://se.mathworks.com/matlabcentral/fileexchange/18801-quadvgk
We here use gaussian kronrod integration already used in gpstuff for evaluating one dimensional integrals.
This is vectorised quadrature which means that several functions can be evaluated at the same time over a grid of
points.
:param f:
:param fmin:
:param fmax:
:param difftol:
:return:
"""
XK = np.array([-0.991455371120813, -0.949107912342759, -0.864864423359769, -0.741531185599394,
-0.586087235467691, -0.405845151377397, -0.207784955007898, 0.,
0.207784955007898, 0.405845151377397, 0.586087235467691,
0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813])
WK = np.array([0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525,
0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728,
0.204432940075298, 0.190350578064785, 0.169004726639267,
0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529])
# 7-point Gaussian weightings
WG = np.array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469,
0.381830050505119, 0.279705391489277, 0.129484966168870])
NK = WK.size
G = np.arange(2,NK,2)
tol1 = 1e-4
tol2 = 1e-4
Subs = np.array([[fmin],[fmax]])
# number of functions to evaluate in the feval vector of functions.
NF = feval(np.zeros(1)).size
Q = np.zeros(NF)
neval = 0
while Subs.size > 0:
Subs = getSubs(Subs,XK)
M = (Subs[1,:] - Subs[0,:]) / 2
C = (Subs[1,:] + Subs[0,:]) / 2
# NM = length(M);
NM = M.size
# x = reshape(XK * M + ones(NK, 1) * C, 1, []);
x = XK[:,None]*M + C
x = x.flatten()
FV = feval(x)
# FV = FV[:,None]
Q1 = np.zeros((NF, NM))
Q2 = np.zeros((NF, NM))
# for n=1:NF
# F = reshape(FV(n,:), NK, []);
# Q1(n,:) = M. * sum((WK * ones(1, NM)). * F);
# Q2(n,:) = M. * sum((WG * ones(1, NM)). * F(G,:));
# end
# for i in range(NF):
# F = FV
# F = F.reshape((NK,-1))
# temp_mat = np.sum(np.multiply(WK[:,None]*np.ones((1,NM)), F),axis=0)
# Q1[i,:] = np.multiply(M, temp_mat)
# temp_mat = np.sum(np.multiply(WG[:,None]*np.ones((1, NM)), F[G-1,:]), axis=0)
# Q2[i,:] = np.multiply(M, temp_mat)
# ind = np.where(np.logical_or(np.max(np.abs(Q1 -Q2) / Q1) < tol1, (Subs[1,:] - Subs[0,:]) <= tol2) > 0)[0]
# Q = Q + np.sum(Q1[:,ind], axis=1)
# np.delete(Subs, ind,axis=1)
Q1 = np.dot(FV.reshape(NF, NK, NM).swapaxes(2,1),WK)*M
Q2 = np.dot(FV.reshape(NF, NK, NM).swapaxes(2,1)[:,:,1::2],WG)*M
#ind = np.nonzero(np.logical_or(np.max(np.abs((Q1-Q2)/Q1), 0) < difftol , M < xtol))[0]
ind = np.nonzero(np.logical_or(np.max(np.abs((Q1-Q2)), 0) < tol1 , (Subs[1,:] - Subs[0,:]) < tol2))[0]
Q = Q + np.sum(Q1[:,ind], axis=1)
Subs = np.delete(Subs, ind, axis=1)
return Q
def quadgk_int(f, fmin=-np.inf, fmax=np.inf, difftol=0.1):
"""
Integrate f from fmin to fmax,
do integration by substitution
x = r / (1-r**2)
when r goes from -1 to 1 , x goes from -inf to inf.
the interval for quadgk function is from -1 to +1, so we transform the space from (-inf,inf) to (-1,1)
:param f:
:param fmin:
:param fmax:
:param difftol:
:return:
"""
difftol = 1e-4
def trans_func(r):
r2 = np.square(r)
x = r / (1-r2)
dx_dr = (1 + r2)/(1-r2)**2
return f(x)*dx_dr
integrand = quadvgk(trans_func, -1., 1., difftol, difftol)
return integrand

View file

@ -16,6 +16,9 @@ recursive-include GPy *.c
recursive-include GPy *.h
recursive-include GPy *.pyx
# LICENSE
include LICENSE.txt
# Testing
#include GPy/testing/baseline/*.png
#include GPy/testing/pickle_test.pickle

View file

@ -76,7 +76,7 @@ If that is the case, it is best to clean the repo and reinstall.
[<img src="https://upload.wikimedia.org/wikipedia/commons/8/8e/OS_X-Logo.svg" height=40px>](http://www.apple.com/osx/)
[<img src="https://upload.wikimedia.org/wikipedia/commons/3/35/Tux.svg" height=40px>](https://en.wikipedia.org/wiki/List_of_Linux_distributions)
Python 2.7, 3.4 and higher
Python 2.7, 3.5 and higher
## Citation

View file

@ -3,12 +3,14 @@ environment:
secure: 8/ZjXFwtd1S7ixd7PJOpptupKKEDhm2da/q3unabJ00=
COVERALLS_REPO_TOKEN:
secure: d3Luic/ESkGaWnZrvWZTKrzO+xaVwJWaRCEP0F+K/9DQGPSRZsJ/Du5g3s4XF+tS
gpy_version: 1.6.2
gpy_version: 1.7.7
matrix:
- PYTHON_VERSION: 2.7
MINICONDA: C:\Miniconda-x64
- PYTHON_VERSION: 3.5
MINICONDA: C:\Miniconda35-x64
- PYTHON_VERSION: 3.6
MINICONDA: C:\Miniconda36-x64
#configuration:
# - Debug
@ -62,21 +64,21 @@ deploy_script:
- echo test >> %USERPROFILE%\\.pypirc
- echo[
- echo [pypi] >> %USERPROFILE%\\.pypirc
- echo username:maxz >> %USERPROFILE%\\.pypirc
- echo password:%pip_access% >> %USERPROFILE%\\.pypirc
- echo username = maxz >> %USERPROFILE%\\.pypirc
- echo password = %pip_access% >> %USERPROFILE%\\.pypirc
- echo[
- echo [test] >> %USERPROFILE%\\.pypirc
- echo repository:https://test.pypi.org/legacy/ >> %USERPROFILE%\\.pypirc
- echo username:maxz >> %USERPROFILE%\\.pypirc
- echo password:%pip_access% >> %USERPROFILE%\\.pypirc
- echo repository = https://testpypi.python.org/pypi >> %USERPROFILE%\\.pypirc
- echo username = maxz >> %USERPROFILE%\\.pypirc
- echo password = %pip_access% >> %USERPROFILE%\\.pypirc
- ps: >-
if ($env:APPVEYOR_REPO_BRANCH -eq 'devel') {
twine upload --skip-existing -r test dist/*
If ($env:APPVEYOR_REPO_BRANCH -eq 'devel') {
echo not deploying on devel # twine upload --skip-existing -r test dist/*
}
elseif ($env:APPVEYOR_REPO_BRANCH -eq 'deploy') {
ElseIf ($env:APPVEYOR_REPO_BRANCH -eq 'deploy') {
twine upload --skip-existing dist/*
}
else {
Else {
echo not deploying on other branches
}

View file

@ -1,5 +1,5 @@
[bumpversion]
current_version = 1.6.2
current_version = 1.7.7
tag = True
commit = True

View file

@ -150,7 +150,7 @@ setup(name = 'GPy',
py_modules = ['GPy.__init__'],
test_suite = 'GPy.testing',
setup_requires = ['numpy>=1.7'],
install_requires = ['numpy>=1.7', 'scipy>=0.16', 'six', 'paramz>=0.6.9'],
install_requires = ['numpy>=1.7', 'scipy>=0.16', 'six', 'paramz>=0.7.4'],
extras_require = {'docs':['sphinx'],
'optional':['mpi4py',
'ipython>=4.0.0',
@ -169,8 +169,8 @@ setup(name = 'GPy',
'Operating System :: Microsoft :: Windows',
'Operating System :: POSIX :: Linux',
'Programming Language :: Python :: 2.7',
'Programming Language :: Python :: 3.3',
'Programming Language :: Python :: 3.5',
'Programming Language :: Python :: 3.6',
'Framework :: IPython',
'Intended Audience :: Science/Research',
'Intended Audience :: Developers',