From 11988eea060404096011b5556dcc89165af05f4a Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Mon, 1 Aug 2022 16:33:47 -0700 Subject: [PATCH 01/12] initial commit --- pymc_experimental/gp/hsgp.py | 110 +++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 pymc_experimental/gp/hsgp.py diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py new file mode 100644 index 000000000..718a6bd5d --- /dev/null +++ b/pymc_experimental/gp/hsgp.py @@ -0,0 +1,110 @@ +import numpy as np +import aesara.tensor as at +import pymc as pm + +from pymc.gp.util import ( + JITTER_DEFAULT, + cholesky, + conditioned_vars, + replace_with_values, +) + + +class ExpQuad(pm.gp.cov.ExpQuad): + def psd(self, omega, ls): + D = len(self.active_dims) + ls = at.ones(D) * ls + c = at.power(at.sqrt(2.0 * np.pi), D) * at.prod(ls) + return c * at.exp(-0.5 * at.dot(omega, ls)) + + +class HSGP(pm.gp.Latent): + def __init__(self, M, c=3/2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0)): + self.M = M + self.c = c + self.L = L + self.D = len(cov_func.active_dims) + super().__init__(mean_func=mean_func, cov_func=cov_func) + + def _evaluate_spd(self, cov_func, omega, L, M, Xsd): + #L = at.as_tensor_variable(L) + #omega = np.pi * at.arange(1, M + 1)[:, None] * at.ones_like(L) / (2.0 * L) + #self.omega = omega + cov, scale = cov_func.factor_list + return scale * cov.psd(omega, cov.ls / Xsd) + + def _standardize(self, X): + self.Xmu = at.mean(X[:, self.cov_func.active_dims], axis=0) + self.Xsd = at.std(X[:, self.cov_func.active_dims], axis=0) + Xz = (X[:, self.cov_func.active_dims] - self.Xmu) / self.Xsd + return Xz + + @staticmethod + def _construct_basis(X, L, D, M): + + S = np.meshgrid(*[np.arange(1, 1 + M) for _ in range(D)]) + S = np.vstack([s.flatten() for s in S]).T + eigvals = at.square((np.pi * S) / (2 * L)) + + m_star = at.power(M, D) + phi_shape = (X.shape[0], m_star) + phi = at.ones(phi_shape) + for d in range(D): + c = 1.0 / np.sqrt(L[d]) + phi *= c * at.sin(at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], m_star) + L[d])) + return eigvals, phi + + def _build_prior(self, name, X, **kwargs): + X = at.as_tensor_variable(X) + + if self.L is None: + Xz = self._standardize(X) + La = at.abs(at.min(Xz, axis=0)) + Lb = at.abs(at.max(Xz, axis=0)) + L = self.c * at.max(at.stack((La, Lb)), axis=0) + self.L = L + + else: + if np.isscalar(self.L): + L = [self.L] + + if len(L) != self.D: + raise ValueError("Must provide one L for each active dimension.") + + Xz = X + self.Xmu, self.Xsd = 0.0, 1.0 + L = at.as_tensor_variable(L) + + eigvals, phi = self._construct_basis(Xz, L, self.D, self.M) + + omega = at.sqrt(eigvals) + self.omega = omega + + spd = self._evaluate_spd(self.cov_func, omega, L, self.M, self.Xsd) + + m_star = at.power(self.M, self.D) + beta = pm.Normal(f'{name}_coeffs_', size=m_star) + f = pm.Deterministic(name, self.mean_func(X) + at.squeeze(at.dot(phi * beta, spd))) + return eigvals, phi, spd, f, beta + + def prior(self, name, X, **kwargs): + eigvals, phi, spd, f, beta = self._build_prior(name, X, **kwargs) + + self.eigvals = eigvals + self.phi = phi + self.spd = spd + + self.f, self.beta = f, beta + return f + + def _build_conditional(self, name, Xnew): + Xznew = (Xnew[:, self.cov_func.active_dims] - self.Xmu) / self.Xsd + eigvals, phi = self._construct_basis(Xznew, self.L, self.D, self.M) + omega = at.sqrt(eigvals) + spd = self._evaluate_spd(self.cov_func, omega, self.L, self.M, self.Xsd) + return self.mean_func(Xnew) + at.squeeze(at.dot(phi * self.beta, spd)) + + def conditional(self, name, Xnew): + # TODO, error if Xnew outside bounds given by L + fnew = self._build_conditional(name, Xnew) + return pm.Deterministic(name, fnew) From 578566c5be8531a4fb9ebab3549445d9e1bf7613 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Mon, 1 Aug 2022 19:39:17 -0700 Subject: [PATCH 02/12] remove hsgp from latent approx file --- pymc_experimental/gp/__init__.py | 3 +- pymc_experimental/gp/latent_approx.py | 102 -------------------------- 2 files changed, 1 insertion(+), 104 deletions(-) diff --git a/pymc_experimental/gp/__init__.py b/pymc_experimental/gp/__init__.py index 133de4dc6..2daf8bc19 100644 --- a/pymc_experimental/gp/__init__.py +++ b/pymc_experimental/gp/__init__.py @@ -1,5 +1,4 @@ from pymc_experimental.gp.latent_approx import ( - HSGP, KarhunenLoeveExpansion, ProjectedProcess, -) +) \ No newline at end of file diff --git a/pymc_experimental/gp/latent_approx.py b/pymc_experimental/gp/latent_approx.py index b0591085a..243508dcd 100644 --- a/pymc_experimental/gp/latent_approx.py +++ b/pymc_experimental/gp/latent_approx.py @@ -60,108 +60,6 @@ def conditional(self, name, Xnew, jitter=1e-6, **kwargs): return pm.MvNormal(name, mu=mu, chol=chol) -class HSGP(pm.gp.Latent): - ## inputs: M, c - - def __init__( - self, n_basis, c=3 / 2, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0) - ): - ## TODO: specify either c or L - self.M = n_basis - self.c = c - super().__init__(mean_func=mean_func, cov_func=cov_func) - - def _validate_cov_func(self, cov_func): - ## TODO: actually validate it. Right now this fails unless cov func is exactly - # in the form eta**2 * pm.gp.cov.Matern12(...) and will error otherwise. - cov, scaling_factor = cov_func.factor_list - return scaling_factor, cov.ls, cov.spectral_density - - def prior(self, name, X, **kwargs): - f, Phi, L, spd, beta, Xmu, Xsd = self._build_prior(name, X, **kwargs) - self.X, self.f = X, f - self.Phi, self.L, self.spd, self.beta = Phi, L, spd, beta - self.Xmu, self.Xsd = Xmu, Xsd - return f - - def _generate_basis(self, X, L): - indices = at.arange(1, self.M + 1) - m1 = (np.pi / (2.0 * L)) * at.tile(L + X, self.M) - m2 = at.diag(indices) - Phi = at.sin(m1 @ m2) / at.sqrt(L) - omega = (np.pi * indices) / (2.0 * L) - return Phi, omega - - def _build_prior(self, name, X, **kwargs): - n_obs = np.shape(X)[0] - - # standardize input scale - X = at.as_tensor_variable(X) - Xmu = at.mean(X, axis=0) - Xsd = at.std(X, axis=0) - Xz = (X - Xmu) / Xsd - - # define L using Xz and c - La = at.abs(at.min(Xz)) # .eval()? - Lb = at.max(Xz) - L = self.c * at.max([La, Lb]) - - # make basis and omega, spectral density - Phi, omega = self._generate_basis(Xz, L) - scale, ls, spectral_density = self._validate_cov_func(self.cov_func) - spd = scale * spectral_density(omega, ls / Xsd).flatten() - - beta = pm.Normal(f"{name}_coeffs_", size=self.M) - f = pm.Deterministic(name, self.mean_func(X) + at.dot(Phi * at.sqrt(spd), beta)) - return f, Phi, L, spd, beta, Xmu, Xsd - - def _build_conditional(self, Xnew, Xmu, Xsd, L, beta): - Xnewz = (Xnew - Xmu) / Xsd - Phi, omega = self._generate_basis(Xnewz, L) - scale, ls, spectral_density = self._validate_cov_func(self.cov_func) - spd = scale * spectral_density(omega, ls / Xsd).flatten() - return self.mean_func(Xnew) + at.dot(Phi * at.sqrt(spd), beta) - - def conditional(self, name, Xnew): - # warn about extrapolation - fnew = self._build_conditional(Xnew, self.Xmu, self.Xsd, self.L, self.beta) - return pm.Deterministic(name, fnew) - - -class ExpQuad(pm.gp.cov.ExpQuad): - @staticmethod - def spectral_density(omega, ls): - # univariate spectral denisty, implement multi - return at.sqrt(2 * np.pi) * ls * at.exp(-0.5 * ls**2 * omega**2) - - -class Matern52(pm.gp.cov.Matern52): - @staticmethod - def spectral_density(omega, ls): - # univariate spectral denisty, implement multi - # https://arxiv.org/pdf/1611.06740.pdf - lam = at.sqrt(5) * (1.0 / ls) - return (16.0 / 3.0) * lam**5 * (1.0 / (lam**2 + omega**2) ** 3) - - -class Matern32(pm.gp.cov.Matern32): - @staticmethod - def spectral_density(omega, ls): - # univariate spectral denisty, implement multi - # https://arxiv.org/pdf/1611.06740.pdf - lam = np.sqrt(3.0) * (1.0 / ls) - return 4.0 * lam**3 * (1.0 / at.square(lam**2 + omega**2)) - - -class Matern12(pm.gp.cov.Matern12): - @staticmethod - def spectral_density(omega, ls): - # univariate spectral denisty, implement multi - # https://arxiv.org/pdf/1611.06740.pdf - lam = 1.0 / ls - return 2.0 * lam * (1.0 / (lam**2 + omega**2)) - - class KarhunenLoeveExpansion(pm.gp.Latent): def __init__( self, From 8cd588895776d20f5eeca7ebc18ff19e1c7800a8 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Thu, 4 Aug 2022 17:25:24 -0700 Subject: [PATCH 03/12] fix init --- pymc_experimental/gp/__init__.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pymc_experimental/gp/__init__.py b/pymc_experimental/gp/__init__.py index 2daf8bc19..4e5499e99 100644 --- a/pymc_experimental/gp/__init__.py +++ b/pymc_experimental/gp/__init__.py @@ -1,4 +1,9 @@ from pymc_experimental.gp.latent_approx import ( KarhunenLoeveExpansion, ProjectedProcess, +) +from pymc_experimental.hsgp import ( + HSGP, + ExpQuad, + Matern52, ) \ No newline at end of file From fb254b8efbf52e0a616b73e090fd76ac8003de87 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Thu, 4 Aug 2022 17:25:40 -0700 Subject: [PATCH 04/12] add matern52 --- pymc_experimental/gp/hsgp.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index 718a6bd5d..277776021 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -15,10 +15,20 @@ def psd(self, omega, ls): D = len(self.active_dims) ls = at.ones(D) * ls c = at.power(at.sqrt(2.0 * np.pi), D) * at.prod(ls) - return c * at.exp(-0.5 * at.dot(omega, ls)) + return c * at.exp(-0.5 * at.dot(omega, ls)) -class HSGP(pm.gp.Latent): +class Matern52(pm.gp.cov.Matern52): + def psd(self, omega, ls): + D = len(self.active_dims) + ls = at.ones(D) * ls + D52 = (D + 5) / 2 + num = at.power(2, D) * at.power(np.pi, D / 2) * at.gamma(D52) * at.power(5, 5 / 2) + den = 0.75 * at.sqrt(np.pi) * at.power(ls, 5) + return (num / den) * at.power(5.0 + at.dot(omega, ls), D52) + + +class HSGP(pm.gp.Base): def __init__(self, M, c=3/2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0)): self.M = M self.c = c @@ -27,9 +37,6 @@ def __init__(self, M, c=3/2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm super().__init__(mean_func=mean_func, cov_func=cov_func) def _evaluate_spd(self, cov_func, omega, L, M, Xsd): - #L = at.as_tensor_variable(L) - #omega = np.pi * at.arange(1, M + 1)[:, None] * at.ones_like(L) / (2.0 * L) - #self.omega = omega cov, scale = cov_func.factor_list return scale * cov.psd(omega, cov.ls / Xsd) @@ -41,7 +48,7 @@ def _standardize(self, X): @staticmethod def _construct_basis(X, L, D, M): - + # S = np.meshgrid(*[np.arange(1, 1 + M) for _ in range(D)]) S = np.vstack([s.flatten() for s in S]).T eigvals = at.square((np.pi * S) / (2 * L)) From c03325a385a696b414bfea5b32b7e421e26ec5c9 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Thu, 4 Aug 2022 20:17:32 -0700 Subject: [PATCH 05/12] get compound kernels working --- pymc_experimental/gp/__init__.py | 2 +- pymc_experimental/gp/hsgp.py | 129 ++++++++++++++++++------------- 2 files changed, 78 insertions(+), 53 deletions(-) diff --git a/pymc_experimental/gp/__init__.py b/pymc_experimental/gp/__init__.py index 4e5499e99..b6b8ff078 100644 --- a/pymc_experimental/gp/__init__.py +++ b/pymc_experimental/gp/__init__.py @@ -2,7 +2,7 @@ KarhunenLoeveExpansion, ProjectedProcess, ) -from pymc_experimental.hsgp import ( +from pymc_experimental.gp.hsgp import ( HSGP, ExpQuad, Matern52, diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index 277776021..d9640f5ea 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -11,105 +11,130 @@ class ExpQuad(pm.gp.cov.ExpQuad): - def psd(self, omega, ls): + def psd(self, omega): D = len(self.active_dims) - ls = at.ones(D) * ls + ls = at.ones(D) * self.ls c = at.power(at.sqrt(2.0 * np.pi), D) * at.prod(ls) return c * at.exp(-0.5 * at.dot(omega, ls)) class Matern52(pm.gp.cov.Matern52): - def psd(self, omega, ls): + def psd(self, omega): D = len(self.active_dims) - ls = at.ones(D) * ls + ls = at.ones(D) * self.ls D52 = (D + 5) / 2 num = at.power(2, D) * at.power(np.pi, D / 2) * at.gamma(D52) * at.power(5, 5 / 2) den = 0.75 * at.sqrt(np.pi) * at.power(ls, 5) return (num / den) * at.power(5.0 + at.dot(omega, ls), D52) -class HSGP(pm.gp.Base): +class HSGP(pm.gp.gp.Base): def __init__(self, M, c=3/2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0)): self.M = M self.c = c self.L = L self.D = len(cov_func.active_dims) + self.m_star = at.power(self.M, self.D) super().__init__(mean_func=mean_func, cov_func=cov_func) - def _evaluate_spd(self, cov_func, omega, L, M, Xsd): - cov, scale = cov_func.factor_list - return scale * cov.psd(omega, cov.ls / Xsd) + def _evaluate_psd(self, omega): + # Unpack and validate covariance function + #cov, scale = self.cov_func.factor_list + + return self.cov_func.psd(omega) + + @staticmethod def _standardize(self, X): - self.Xmu = at.mean(X[:, self.cov_func.active_dims], axis=0) - self.Xsd = at.std(X[:, self.cov_func.active_dims], axis=0) - Xz = (X[:, self.cov_func.active_dims] - self.Xmu) / self.Xsd - return Xz + Xmu = at.mean(X[:, self.cov_func.active_dims], axis=0) + Xsd = at.std(X[:, self.cov_func.active_dims], axis=0) + Xz = (X[:, self.cov_func.active_dims] - Xmu) / Xsd + return Xz, Xmu, Xsd - @staticmethod - def _construct_basis(X, L, D, M): - # - S = np.meshgrid(*[np.arange(1, 1 + M) for _ in range(D)]) + def _construct_basis(self, X): + """ Construct the set of basis vectors and associated eigenvalue array. + """ + S = np.meshgrid(*[np.arange(1, 1 + self.M) for _ in range(self.D)]) S = np.vstack([s.flatten() for s in S]).T - eigvals = at.square((np.pi * S) / (2 * L)) + eigvals = at.square((np.pi * S) / (2 * self.L)) - m_star = at.power(M, D) - phi_shape = (X.shape[0], m_star) + phi_shape = (X.shape[0], self.m_star) phi = at.ones(phi_shape) - for d in range(D): - c = 1.0 / np.sqrt(L[d]) - phi *= c * at.sin(at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], m_star) + L[d])) + for d in range(self.D): + c = 1.0 / np.sqrt(self.L[d]) + phi *= c * at.sin(at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], self.m_star) + self.L[d])) return eigvals, phi def _build_prior(self, name, X, **kwargs): X = at.as_tensor_variable(X) if self.L is None: - Xz = self._standardize(X) - La = at.abs(at.min(Xz, axis=0)) - Lb = at.abs(at.max(Xz, axis=0)) - L = self.c * at.max(at.stack((La, Lb)), axis=0) - self.L = L + + # Define new L based on c and X range + La = at.abs(at.min(X, axis=0)) + Lb = at.abs(at.max(X, axis=0)) + self.L = self.c * at.max(at.stack((La, Lb)), axis=0) else: + # If L is passed as a scalar, put it into a one-element list. if np.isscalar(self.L): - L = [self.L] - - if len(L) != self.D: - raise ValueError("Must provide one L for each active dimension.") - - Xz = X - self.Xmu, self.Xsd = 0.0, 1.0 - L = at.as_tensor_variable(L) - - eigvals, phi = self._construct_basis(Xz, L, self.D, self.M) - + self.L = [self.L] + + if len(self.L) != self.D: + raise ValueError( + ( + "Must provide one L for each active dimension. `len(L)` must ", + "equal the number of active dimensions of your covariance function." + ) + ) + + ## If L is provided, don't rescale X + self.L = at.as_tensor_variable(self.L) + + # Construct basis and eigenvalues + eigvals, phi = self._construct_basis(X) omega = at.sqrt(eigvals) - self.omega = omega - - spd = self._evaluate_spd(self.cov_func, omega, L, self.M, self.Xsd) - - m_star = at.power(self.M, self.D) - beta = pm.Normal(f'{name}_coeffs_', size=m_star) - f = pm.Deterministic(name, self.mean_func(X) + at.squeeze(at.dot(phi * beta, spd))) - return eigvals, phi, spd, f, beta + psd = self._evaluate_psd(omega) + beta = pm.Normal(f'{name}_coeffs_', size=self.m_star) + f = pm.Deterministic(name, self.mean_func(X) + at.squeeze(at.dot(phi * beta, psd))) + return eigvals, phi, psd, f, beta + @property + def basis(self): + try: + return self.phi + except AttributeError: + raise RuntimeError("Must construct the prior first by calling `.prior.") + + @property + def power_spectral_density(self): + try: + return self.psd + except AttributeError: + raise RuntimeError("Must construct the prior first by calling `.prior.") + + @property + def omega(self): + try: + return at.sqrt(self.eigvals) + except AttributeError: + raise RuntimeError("Must construct the prior first by calling `.prior.") + def prior(self, name, X, **kwargs): - eigvals, phi, spd, f, beta = self._build_prior(name, X, **kwargs) + eigvals, phi, psd, f, beta = self._build_prior(name, X, **kwargs) self.eigvals = eigvals self.phi = phi - self.spd = spd + self.psd = psd self.f, self.beta = f, beta return f def _build_conditional(self, name, Xnew): - Xznew = (Xnew[:, self.cov_func.active_dims] - self.Xmu) / self.Xsd - eigvals, phi = self._construct_basis(Xznew, self.L, self.D, self.M) + eigvals, phi = self._construct_basis(Xnew) omega = at.sqrt(eigvals) - spd = self._evaluate_spd(self.cov_func, omega, self.L, self.M, self.Xsd) - return self.mean_func(Xnew) + at.squeeze(at.dot(phi * self.beta, spd)) + psd = self._evaluate_psd(omega) + return self.mean_func(Xnew) + at.squeeze(at.dot(phi * self.beta, psd)) def conditional(self, name, Xnew): # TODO, error if Xnew outside bounds given by L From 0d48e357699843acbc373309b42f54c3f6b22c91 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Thu, 4 Aug 2022 23:54:42 -0700 Subject: [PATCH 06/12] tidy a bit --- pymc_experimental/gp/hsgp.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index d9640f5ea..1bd68d1a7 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -121,14 +121,8 @@ def omega(self): raise RuntimeError("Must construct the prior first by calling `.prior.") def prior(self, name, X, **kwargs): - eigvals, phi, psd, f, beta = self._build_prior(name, X, **kwargs) - - self.eigvals = eigvals - self.phi = phi - self.psd = psd - - self.f, self.beta = f, beta - return f + self.eigvals, self.phi, self.psd, self.f, self.beta = self._build_prior(name, X, **kwargs) + return self.f def _build_conditional(self, name, Xnew): eigvals, phi = self._construct_basis(Xnew) @@ -137,6 +131,5 @@ def _build_conditional(self, name, Xnew): return self.mean_func(Xnew) + at.squeeze(at.dot(phi * self.beta, psd)) def conditional(self, name, Xnew): - # TODO, error if Xnew outside bounds given by L fnew = self._build_conditional(name, Xnew) - return pm.Deterministic(name, fnew) + return pm.Deterministic(name, fnew) \ No newline at end of file From 11fd5ed049b1e295ee75d373897862f8ca8bea5d Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Sun, 7 Aug 2022 11:21:58 -0700 Subject: [PATCH 07/12] clean up, move cov function psds to pymc pr --- pymc_experimental/gp/hsgp.py | 52 +++++++----------------------------- 1 file changed, 10 insertions(+), 42 deletions(-) diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index 1bd68d1a7..5adcae252 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -8,59 +8,26 @@ conditioned_vars, replace_with_values, ) - - -class ExpQuad(pm.gp.cov.ExpQuad): - def psd(self, omega): - D = len(self.active_dims) - ls = at.ones(D) * self.ls - c = at.power(at.sqrt(2.0 * np.pi), D) * at.prod(ls) - return c * at.exp(-0.5 * at.dot(omega, ls)) - - -class Matern52(pm.gp.cov.Matern52): - def psd(self, omega): - D = len(self.active_dims) - ls = at.ones(D) * self.ls - D52 = (D + 5) / 2 - num = at.power(2, D) * at.power(np.pi, D / 2) * at.gamma(D52) * at.power(5, 5 / 2) - den = 0.75 * at.sqrt(np.pi) * at.power(ls, 5) - return (num / den) * at.power(5.0 + at.dot(omega, ls), D52) class HSGP(pm.gp.gp.Base): def __init__(self, M, c=3/2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0)): + super().__init__(mean_func=mean_func, cov_func=cov_func) self.M = M self.c = c self.L = L - self.D = len(cov_func.active_dims) - self.m_star = at.power(self.M, self.D) - super().__init__(mean_func=mean_func, cov_func=cov_func) - - def _evaluate_psd(self, omega): - - # Unpack and validate covariance function - #cov, scale = self.cov_func.factor_list + self.m_star = at.power(self.M, self.cov_func.D) - return self.cov_func.psd(omega) - - @staticmethod - def _standardize(self, X): - Xmu = at.mean(X[:, self.cov_func.active_dims], axis=0) - Xsd = at.std(X[:, self.cov_func.active_dims], axis=0) - Xz = (X[:, self.cov_func.active_dims] - Xmu) / Xsd - return Xz, Xmu, Xsd - def _construct_basis(self, X): """ Construct the set of basis vectors and associated eigenvalue array. """ - S = np.meshgrid(*[np.arange(1, 1 + self.M) for _ in range(self.D)]) + S = np.meshgrid(*[np.arange(1, 1 + self.M) for _ in range(self.cov_func.D)]) S = np.vstack([s.flatten() for s in S]).T eigvals = at.square((np.pi * S) / (2 * self.L)) phi_shape = (X.shape[0], self.m_star) phi = at.ones(phi_shape) - for d in range(self.D): + for d in range(self.cov_func.D): c = 1.0 / np.sqrt(self.L[d]) phi *= c * at.sin(at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], self.m_star) + self.L[d])) return eigvals, phi @@ -80,7 +47,8 @@ def _build_prior(self, name, X, **kwargs): if np.isscalar(self.L): self.L = [self.L] - if len(self.L) != self.D: + # Make sure L has the right dimension + if len(self.L) != self.cov_func.D: raise ValueError( ( "Must provide one L for each active dimension. `len(L)` must ", @@ -94,9 +62,9 @@ def _build_prior(self, name, X, **kwargs): # Construct basis and eigenvalues eigvals, phi = self._construct_basis(X) omega = at.sqrt(eigvals) - psd = self._evaluate_psd(omega) + psd = self.cov_func.psd(omega) beta = pm.Normal(f'{name}_coeffs_', size=self.m_star) - f = pm.Deterministic(name, self.mean_func(X) + at.squeeze(at.dot(phi * beta, psd))) + f = pm.Deterministic(name, self.mean_func(X) + at.squeeze(at.dot(phi, beta * psd))) return eigvals, phi, psd, f, beta @property @@ -127,8 +95,8 @@ def prior(self, name, X, **kwargs): def _build_conditional(self, name, Xnew): eigvals, phi = self._construct_basis(Xnew) omega = at.sqrt(eigvals) - psd = self._evaluate_psd(omega) - return self.mean_func(Xnew) + at.squeeze(at.dot(phi * self.beta, psd)) + psd = self.cov_func.psd(omega) + return self.mean_func(Xnew) + at.squeeze(at.dot(phi, self.beta * psd)) def conditional(self, name, Xnew): fnew = self._build_conditional(name, Xnew) From 8323f189020a1ff826fdc265a2a26dc6de2410b5 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Sun, 7 Aug 2022 11:22:10 -0700 Subject: [PATCH 08/12] fix init --- pymc_experimental/gp/__init__.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/pymc_experimental/gp/__init__.py b/pymc_experimental/gp/__init__.py index b6b8ff078..f3d5f8869 100644 --- a/pymc_experimental/gp/__init__.py +++ b/pymc_experimental/gp/__init__.py @@ -2,8 +2,4 @@ KarhunenLoeveExpansion, ProjectedProcess, ) -from pymc_experimental.gp.hsgp import ( - HSGP, - ExpQuad, - Matern52, -) \ No newline at end of file +from pymc_experimental.gp.hsgp import HSGP \ No newline at end of file From b866404636f02938447b92222e876fd2a7f9fbb3 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Sun, 7 Aug 2022 11:23:12 -0700 Subject: [PATCH 09/12] run precommit --- pymc_experimental/gp/__init__.py | 7 ++-- pymc_experimental/gp/hsgp.py | 60 +++++++++++++++----------------- 2 files changed, 30 insertions(+), 37 deletions(-) diff --git a/pymc_experimental/gp/__init__.py b/pymc_experimental/gp/__init__.py index f3d5f8869..e9f043778 100644 --- a/pymc_experimental/gp/__init__.py +++ b/pymc_experimental/gp/__init__.py @@ -1,5 +1,2 @@ -from pymc_experimental.gp.latent_approx import ( - KarhunenLoeveExpansion, - ProjectedProcess, -) -from pymc_experimental.gp.hsgp import HSGP \ No newline at end of file +from pymc_experimental.gp.hsgp import HSGP +from pymc_experimental.gp.latent_approx import KarhunenLoeveExpansion, ProjectedProcess diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index 5adcae252..28a21a048 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -1,26 +1,20 @@ -import numpy as np import aesara.tensor as at +import numpy as np import pymc as pm -from pymc.gp.util import ( - JITTER_DEFAULT, - cholesky, - conditioned_vars, - replace_with_values, -) - class HSGP(pm.gp.gp.Base): - def __init__(self, M, c=3/2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0)): + def __init__( + self, M, c=3 / 2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0) + ): super().__init__(mean_func=mean_func, cov_func=cov_func) self.M = M self.c = c self.L = L self.m_star = at.power(self.M, self.cov_func.D) - - def _construct_basis(self, X): - """ Construct the set of basis vectors and associated eigenvalue array. - """ + + def _construct_basis(self, X): + """Construct the set of basis vectors and associated eigenvalue array.""" S = np.meshgrid(*[np.arange(1, 1 + self.M) for _ in range(self.cov_func.D)]) S = np.vstack([s.flatten() for s in S]).T eigvals = at.square((np.pi * S) / (2 * self.L)) @@ -29,75 +23,77 @@ def _construct_basis(self, X): phi = at.ones(phi_shape) for d in range(self.cov_func.D): c = 1.0 / np.sqrt(self.L[d]) - phi *= c * at.sin(at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], self.m_star) + self.L[d])) - return eigvals, phi + phi *= c * at.sin( + at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], self.m_star) + self.L[d]) + ) + return eigvals, phi def _build_prior(self, name, X, **kwargs): X = at.as_tensor_variable(X) - + if self.L is None: - + # Define new L based on c and X range La = at.abs(at.min(X, axis=0)) Lb = at.abs(at.max(X, axis=0)) self.L = self.c * at.max(at.stack((La, Lb)), axis=0) - + else: # If L is passed as a scalar, put it into a one-element list. if np.isscalar(self.L): self.L = [self.L] - + # Make sure L has the right dimension if len(self.L) != self.cov_func.D: raise ValueError( ( "Must provide one L for each active dimension. `len(L)` must ", - "equal the number of active dimensions of your covariance function." + "equal the number of active dimensions of your covariance function.", ) ) - + ## If L is provided, don't rescale X self.L = at.as_tensor_variable(self.L) - + # Construct basis and eigenvalues eigvals, phi = self._construct_basis(X) omega = at.sqrt(eigvals) psd = self.cov_func.psd(omega) - beta = pm.Normal(f'{name}_coeffs_', size=self.m_star) + beta = pm.Normal(f"{name}_coeffs_", size=self.m_star) f = pm.Deterministic(name, self.mean_func(X) + at.squeeze(at.dot(phi, beta * psd))) return eigvals, phi, psd, f, beta - + @property def basis(self): try: return self.phi except AttributeError: - raise RuntimeError("Must construct the prior first by calling `.prior.") - + raise RuntimeError("Must construct the prior first by calling `.prior.") + @property def power_spectral_density(self): try: return self.psd except AttributeError: - raise RuntimeError("Must construct the prior first by calling `.prior.") - + raise RuntimeError("Must construct the prior first by calling `.prior.") + @property def omega(self): try: return at.sqrt(self.eigvals) except AttributeError: - raise RuntimeError("Must construct the prior first by calling `.prior.") + raise RuntimeError("Must construct the prior first by calling `.prior.") def prior(self, name, X, **kwargs): self.eigvals, self.phi, self.psd, self.f, self.beta = self._build_prior(name, X, **kwargs) return self.f - + def _build_conditional(self, name, Xnew): eigvals, phi = self._construct_basis(Xnew) omega = at.sqrt(eigvals) psd = self.cov_func.psd(omega) return self.mean_func(Xnew) + at.squeeze(at.dot(phi, self.beta * psd)) - + def conditional(self, name, Xnew): fnew = self._build_conditional(name, Xnew) - return pm.Deterministic(name, fnew) \ No newline at end of file + return pm.Deterministic(name, fnew) From a432ae23822dd7266e073e38451d647afa8ca792 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Fri, 12 Aug 2022 13:34:59 -0700 Subject: [PATCH 10/12] hsgp working --- pymc_experimental/gp/hsgp.py | 121 ++++++++++++++++------------------- 1 file changed, 55 insertions(+), 66 deletions(-) diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index 28a21a048..0bac46794 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -5,92 +5,81 @@ class HSGP(pm.gp.gp.Base): def __init__( - self, M, c=3 / 2, L=None, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0) + self, + n_basis, + c=3 / 2, + L=None, + *, + mean_func=pm.gp.mean.Zero(), + cov_func=pm.gp.cov.Constant(0.0), ): - super().__init__(mean_func=mean_func, cov_func=cov_func) - self.M = M - self.c = c - self.L = L - self.m_star = at.power(self.M, self.cov_func.D) + arg_err_msg = ( + "`n_basis` and L, if provided, must be lists or tuples, with one element per active " + "dimension." + ) + try: + if len(n_basis) != cov_func.D: + raise ValueError(arg_err_msg) + except TypeError as e: + raise ValueError(arg_err_msg) from e - def _construct_basis(self, X): - """Construct the set of basis vectors and associated eigenvalue array.""" - S = np.meshgrid(*[np.arange(1, 1 + self.M) for _ in range(self.cov_func.D)]) - S = np.vstack([s.flatten() for s in S]).T - eigvals = at.square((np.pi * S) / (2 * self.L)) + if L is not None and len(L) != cov_func.D: + raise ValueError(arg_err_msg) + + self.M = n_basis + self.L = L + self.c = c + self.D = cov_func.D - phi_shape = (X.shape[0], self.m_star) - phi = at.ones(phi_shape) - for d in range(self.cov_func.D): - c = 1.0 / np.sqrt(self.L[d]) - phi *= c * at.sin( - at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], self.m_star) + self.L[d]) - ) - return eigvals, phi + super().__init__(mean_func=mean_func, cov_func=cov_func) - def _build_prior(self, name, X, **kwargs): - X = at.as_tensor_variable(X) + def __add__(self, other): + raise NotImplementedError("Additive HSGP's isn't supported ") + def _set_boundary(self, X): if self.L is None: - # Define new L based on c and X range La = at.abs(at.min(X, axis=0)) Lb = at.abs(at.max(X, axis=0)) self.L = self.c * at.max(at.stack((La, Lb)), axis=0) - else: - # If L is passed as a scalar, put it into a one-element list. - if np.isscalar(self.L): - self.L = [self.L] - - # Make sure L has the right dimension - if len(self.L) != self.cov_func.D: - raise ValueError( - ( - "Must provide one L for each active dimension. `len(L)` must ", - "equal the number of active dimensions of your covariance function.", - ) - ) - - ## If L is provided, don't rescale X self.L = at.as_tensor_variable(self.L) - # Construct basis and eigenvalues - eigvals, phi = self._construct_basis(X) + @staticmethod + def _eigendecomposition(X, L, M, D): + """Construct the eigenvalues and eigenfunctions of the Laplace operator.""" + m_star = at.prod(M) + S = np.meshgrid(*[np.arange(1, 1 + M[d]) for d in range(D)]) + S = np.vstack([s.flatten() for s in S]).T + eigvals = at.square((np.pi * S) / (2 * L)) + phi = at.ones((X.shape[0], m_star)) + for d in range(D): + c = 1.0 / np.sqrt(L[d]) + phi *= c * at.sin(at.sqrt(eigvals[:, d]) * (at.tile(X[:, d][:, None], m_star) + L[d])) omega = at.sqrt(eigvals) - psd = self.cov_func.psd(omega) - beta = pm.Normal(f"{name}_coeffs_", size=self.m_star) - f = pm.Deterministic(name, self.mean_func(X) + at.squeeze(at.dot(phi, beta * psd))) - return eigvals, phi, psd, f, beta + return omega, phi, m_star - @property - def basis(self): - try: - return self.phi - except AttributeError: - raise RuntimeError("Must construct the prior first by calling `.prior.") - - @property - def power_spectral_density(self): - try: - return self.psd - except AttributeError: - raise RuntimeError("Must construct the prior first by calling `.prior.") - - @property - def omega(self): - try: - return at.sqrt(self.eigvals) - except AttributeError: - raise RuntimeError("Must construct the prior first by calling `.prior.") + def approx_K(self, X): + # construct eigenvectors after slicing X + X, _ = self.cov_func._slice(X) + omega, phi, _ = self._eigendecomposition(X, self.L, self.M, self.D) + psd = self.cov_func.psd(omega) + return at.dot(phi * psd, at.transpose(phi)) def prior(self, name, X, **kwargs): - self.eigvals, self.phi, self.psd, self.f, self.beta = self._build_prior(name, X, **kwargs) + X, _ = self.cov_func._slice(X) + self._set_boundary(X) + omega, phi, m_star = self._eigendecomposition(X, self.L, self.M, self.D) + psd = self.cov_func.psd(omega) + self.beta = pm.Normal(f"{name}_coeffs_", size=m_star) + self.f = pm.Deterministic( + name, self.mean_func(X) + at.squeeze(at.dot(phi, self.beta * psd)) + ) return self.f def _build_conditional(self, name, Xnew): - eigvals, phi = self._construct_basis(Xnew) - omega = at.sqrt(eigvals) + Xnew, _ = self.cov_func._slice(Xnew) + omega, phi, _ = self._eigendecomposition(Xnew, self.L, self.M, self.D) psd = self.cov_func.psd(omega) return self.mean_func(Xnew) + at.squeeze(at.dot(phi, self.beta * psd)) From 2aa4091a4dcd84cd5ceba7d8d16d4ec93fa85ca9 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Fri, 12 Aug 2022 23:28:55 -0700 Subject: [PATCH 11/12] add docstring --- pymc_experimental/gp/hsgp.py | 70 ++++++++++++++++++++++++++++++++++-- 1 file changed, 68 insertions(+), 2 deletions(-) diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index 0bac46794..8eb46d12f 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -4,11 +4,78 @@ class HSGP(pm.gp.gp.Base): + R""" + Hilbert Space Gaussian process + + The `gp.HSGP` class is an implementation of the Hilbert Space Gaussian process. This + approximation is a linear model that uses a fixed set of basis vectors, whose coeficients are + random functions of a stationary covariance function's power spectral density. Like + `gp.Latent`, it does not assume a Gaussian noise model and can be used with any likelihood or as + a component anywhere within a model. Also like `gp.Latent`, it has `prior` and `conditional` + methods. It additonally has an `approx_K` method which returns the approximate covariance + matrix. It supports a limited subset of additive covariances. + + Parameters + ---------- + n_basis: list + The number of basis vectors to use for each active dimension (covariance parameter + `active_dim`). + L: list + The boundary of the space for each `active_dim`. It is called the boundary condition. + Choose L such that the domain `[-L, L]` contains all points in the column of X given by the + `active_dim`. + c: 3/2 + The proportion extension factor. Used to construct L from X. Defined as `S = max|X|` such + that `X` is in `[-S, S]`. `L` is the calculated as `c * S`. One of `c` or `L` must be + provided. Further information can be found in Ruitort-Mayol et. al. + cov_func: None, 2D array, or instance of Covariance + The covariance function. Defaults to zero. + mean_func: None, instance of Mean + The mean function. Defaults to zero. + + Examples + -------- + .. code:: python + + # A three dimensional column vector of inputs. + X = np.random.randn(100, 3) + + with pm.Model() as model: + # Specify the covariance function. Three input dimensions, but we only want to use the + # last two. + cov_func = pm.gp.cov.ExpQuad(3, ls=0.1, active_dims=[1, 2]) + + # Specify the HSGP. Use 10 basis vectors across each active dimension, [1, 2] for a + # total of 10 * 10 = 100. The input X is normally distributed, so use a boundary + # condition that should easily contain all the points, from -6 to 6 in each dimension. + gp = pmx.gp.HSGP(n_basis=[10, 10], L=[6, 6], cov_func=cov_func) + + # Place a GP prior over the function f. + f = gp.prior("f", X=X) + + ... + + # After fitting or sampling, specify the distribution + # at new points with .conditional + Xnew = np.linspace(-1, 2, 50)[:, None] + + with model: + fcond = gp.conditional("fcond", Xnew=Xnew) + + References + ---------- + - Ruitort-Mayol, G., and Anderson, M., and Solin, A., and Vehtari, A. (2022). Practical + Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming + + - Solin, A., Sarkka, S. (2019) Hilbert Space Methods for Reduced-Rank Gaussian Process + Regression. + """ + def __init__( self, n_basis, - c=3 / 2, L=None, + c=3 / 2, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0), @@ -60,7 +127,6 @@ def _eigendecomposition(X, L, M, D): return omega, phi, m_star def approx_K(self, X): - # construct eigenvectors after slicing X X, _ = self.cov_func._slice(X) omega, phi, _ = self._eigendecomposition(X, self.L, self.M, self.D) psd = self.cov_func.psd(omega) From cd6a8a1c01c5d2668424f6e6613a85db375f9f06 Mon Sep 17 00:00:00 2001 From: Bill Engels Date: Tue, 30 Aug 2022 17:44:39 -0700 Subject: [PATCH 12/12] fix docs, consistent notation with paper --- pymc_experimental/gp/hsgp.py | 79 +++++++++++++++++++++++++++--------- 1 file changed, 60 insertions(+), 19 deletions(-) diff --git a/pymc_experimental/gp/hsgp.py b/pymc_experimental/gp/hsgp.py index 8eb46d12f..e77298bd7 100644 --- a/pymc_experimental/gp/hsgp.py +++ b/pymc_experimental/gp/hsgp.py @@ -1,3 +1,5 @@ +import warnings + import aesara.tensor as at import numpy as np import pymc as pm @@ -15,16 +17,19 @@ class HSGP(pm.gp.gp.Base): methods. It additonally has an `approx_K` method which returns the approximate covariance matrix. It supports a limited subset of additive covariances. + For information on choosing appropriate `m`, `L`, and `c`, refer Ruitort-Mayol et. al. or to the + pymc examples documentation. + Parameters ---------- - n_basis: list + m: list The number of basis vectors to use for each active dimension (covariance parameter `active_dim`). L: list The boundary of the space for each `active_dim`. It is called the boundary condition. Choose L such that the domain `[-L, L]` contains all points in the column of X given by the `active_dim`. - c: 3/2 + c: 1.5 The proportion extension factor. Used to construct L from X. Defined as `S = max|X|` such that `X` is in `[-S, S]`. `L` is the calculated as `c * S`. One of `c` or `L` must be provided. Further information can be found in Ruitort-Mayol et. al. @@ -73,19 +78,19 @@ class HSGP(pm.gp.gp.Base): def __init__( self, - n_basis, + m, L=None, - c=3 / 2, + c=1.5, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0), ): arg_err_msg = ( - "`n_basis` and L, if provided, must be lists or tuples, with one element per active " + "`m` and L, if provided, must be lists or tuples, with one element per active " "dimension." ) try: - if len(n_basis) != cov_func.D: + if len(m) != cov_func.D: raise ValueError(arg_err_msg) except TypeError as e: raise ValueError(arg_err_msg) from e @@ -93,7 +98,13 @@ def __init__( if L is not None and len(L) != cov_func.D: raise ValueError(arg_err_msg) - self.M = n_basis + if L is None and c < 1.2: + warnings.warn( + "Most applications will require a `c >= 1.2` for accuracy at the boundaries of the " + "domain." + ) + + self.m = m self.L = L self.c = c self.D = cov_func.D @@ -101,9 +112,10 @@ def __init__( super().__init__(mean_func=mean_func, cov_func=cov_func) def __add__(self, other): - raise NotImplementedError("Additive HSGP's isn't supported ") + raise NotImplementedError("Additive HSGPs aren't supported ") def _set_boundary(self, X): + """Make L from X and c if L is not passed in.""" if self.L is None: # Define new L based on c and X range La = at.abs(at.min(X, axis=0)) @@ -113,10 +125,10 @@ def _set_boundary(self, X): self.L = at.as_tensor_variable(self.L) @staticmethod - def _eigendecomposition(X, L, M, D): + def _eigendecomposition(X, L, m, D): """Construct the eigenvalues and eigenfunctions of the Laplace operator.""" - m_star = at.prod(M) - S = np.meshgrid(*[np.arange(1, 1 + M[d]) for d in range(D)]) + m_star = at.prod(m) + S = np.meshgrid(*[np.arange(1, 1 + m[d]) for d in range(D)]) S = np.vstack([s.flatten() for s in S]).T eigvals = at.square((np.pi * S) / (2 * L)) phi = at.ones((X.shape[0], m_star)) @@ -126,29 +138,58 @@ def _eigendecomposition(X, L, M, D): omega = at.sqrt(eigvals) return omega, phi, m_star - def approx_K(self, X): + def approx_K(self, X, L, m): + """A helper function which gives the approximate kernel or covariance matrix K. This can be + helpful when trying to see how well an approximation may work. + """ X, _ = self.cov_func._slice(X) - omega, phi, _ = self._eigendecomposition(X, self.L, self.M, self.D) + omega, phi, _ = self._eigendecomposition(X, self.L, self.m, self.cov_func.D) psd = self.cov_func.psd(omega) return at.dot(phi * psd, at.transpose(phi)) - def prior(self, name, X, **kwargs): + def prior(self, name, X, dims=None): + R""" + Returns the (approximate) GP prior distribution evaluated over the input locations `X`. + + Parameters + ---------- + name: string + Name of the random variable + X: array-like + Function input values. + dims: None + Dimension name for the GP random variable. + """ + X, _ = self.cov_func._slice(X) self._set_boundary(X) - omega, phi, m_star = self._eigendecomposition(X, self.L, self.M, self.D) + omega, phi, m_star = self._eigendecomposition(X, self.L, self.m, self.D) psd = self.cov_func.psd(omega) self.beta = pm.Normal(f"{name}_coeffs_", size=m_star) self.f = pm.Deterministic( - name, self.mean_func(X) + at.squeeze(at.dot(phi, self.beta * psd)) + name, self.mean_func(X) + at.squeeze(at.dot(phi, self.beta * psd)), dims ) return self.f def _build_conditional(self, name, Xnew): Xnew, _ = self.cov_func._slice(Xnew) - omega, phi, _ = self._eigendecomposition(Xnew, self.L, self.M, self.D) + omega, phi, _ = self._eigendecomposition(Xnew, self.L, self.m, self.D) psd = self.cov_func.psd(omega) return self.mean_func(Xnew) + at.squeeze(at.dot(phi, self.beta * psd)) - def conditional(self, name, Xnew): + def conditional(self, name, Xnew, dims=None): + R""" + Returns the (approximate) conditional distribution evaluated over new input locations + `Xnew`. + + Parameters + ---------- + name: string + Name of the random variable + Xnew: array-like + Function input values. + dims: None + Dimension name for the GP random variable. + """ fnew = self._build_conditional(name, Xnew) - return pm.Deterministic(name, fnew) + return pm.Deterministic(name, fnew, dims)