Source code for glimix_core.ggp._expfam

import warnings

from liknorm import LikNormMachine
from numpy import ascontiguousarray, sign
from numpy.linalg import LinAlgError
from optimix import Function

from .._ep import EP
from .._util import check_outcome


[docs]class ExpFamGP(Function): r"""Expectation Propagation for Generalised Gaussian Processes. Parameters ---------- y : array_like Outcome variable. lik_name : str Likelihood name. mean : function Mean function. (Refer to :doc:`mean`.) cov : function Covariance function. (Refer to :doc:`cov`.) Example ------- .. doctest:: >>> from numpy.random import RandomState >>> >>> from glimix_core.example import offset_mean >>> from glimix_core.example import linear_eye_cov >>> from glimix_core.ggp import ExpFamGP >>> from glimix_core.lik import BernoulliProdLik >>> from glimix_core.link import LogitLink >>> from glimix_core.random import GGPSampler >>> >>> random = RandomState(1) >>> >>> lik = BernoulliProdLik(LogitLink()) >>> mean = offset_mean() >>> cov = linear_eye_cov() >>> >>> y = GGPSampler(lik, mean, cov).sample(random) >>> >>> ggp = ExpFamGP(y, 'bernoulli', mean, cov) >>> print('Before: %.4f' % ggp.lml()) Before: -6.9802 >>> ggp.fit(verbose=False) >>> print('After: %.2f' % ggp.lml()) After: -6.55 """ def __init__(self, y, lik, mean, cov): if isinstance(y, tuple): n = len(y[0]) else: n = len(y) Function.__init__(self, "ExpFamGP", composite=[mean, cov]) if not isinstance(lik, (tuple, list)): lik = (lik,) self._lik = (lik[0].lower(),) + tuple(ascontiguousarray(i) for i in lik[1:]) self._y = check_outcome(y, self._lik) self._mean = mean self._cov = cov self._ep = EP(n) self._ep.set_compute_moments(self.compute_moments) self._mean_value = None self._cov_value = None self._machine = LikNormMachine(lik[0], 500) def __del__(self): if hasattr(self, "_machine"): self._machine.finish()
[docs] def fit(self, verbose=True, factr=1e5, pgtol=1e-7): r"""Maximise the marginal likelihood. Parameters ---------- verbose : bool ``True`` for progress output; ``False`` otherwise. Defaults to ``True``. factr : float, optional The iteration stops when ``(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps``, where ``eps`` is the machine precision. pgtol : float, optional The iteration will stop when ``max{|proj g_i | i = 1, ..., n} <= pgtol`` where ``pg_i`` is the i-th component of the projected gradient. Notes ----- Please, refer to :func:`scipy.optimize.fmin_l_bfgs_b` for further information about ``factr`` and ``pgtol``. """ self._maximize(verbose=verbose, factr=factr, pgtol=pgtol)
[docs] def lml(self): r"""Log of the marginal likelihood. Returns ------- float :math:`\log p(\mathbf y)` """ return self.value()
def compute_moments(self, eta, tau, moments): y = (self._y,) + self._lik[1:] self._machine.moments(y, eta, tau, moments) def value(self): from numpy_sugar import epsilon from numpy_sugar.linalg import economic_qs mean = self._mean.value() cov = self._cov.value() try: self._ep.set_prior(mean, dict(QS=economic_qs(cov))) lml = self._ep.lml() except (ValueError, LinAlgError) as e: warnings.warn(str(e), RuntimeWarning) lml = -1 / epsilon.small return lml def gradient(self): from numpy_sugar import epsilon from numpy_sugar.linalg import economic_qs mean = self._mean.value() cov = self._cov.value() gmean = self._mean.gradient() gcov = self._cov.gradient() try: self._ep.set_prior(mean, dict(QS=economic_qs(cov))) grad = dict() for n, g in iter(gmean.items()): grad["ExpFamGP[0]." + n] = self._ep.lml_derivative_over_mean(g) for n, g in iter(gcov.items()): QS = economic_qs(g) grad["ExpFamGP[1]." + n] = self._ep.lml_derivative_over_cov(QS) return grad except (ValueError, LinAlgError) as e: warnings.warn(str(e), RuntimeWarning) v = self._variables.select(fixed=False) return {i: -sign(v.get(i).value) / epsilon.small for i in v}