Add drop_first option to HSGPPeriodic
#7115
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -71,18 +71,17 @@ | |
| Xs: TensorLike, | ||
| period: TensorLike, | ||
| m: int, | ||
| ): | ||
| """ | ||
| Calculate basis vectors for the cosine series expansion of the periodic covariance function. | ||
| These are derived from the Taylor series representation of the covariance. | ||
| """ | ||
| w0 = (2 * pt.pi) / period # angular frequency defining the periodicity | ||
| m1 = pt.tile(w0 * Xs, m) | ||
| m2 = pt.diag(pt.arange(0, m, 1)) | ||
| mw0x = m1 @ m2 | ||
| phi_cos = pt.cos(mw0x) | ||
| phi_sin = pt.sin(mw0x) | ||
| return phi_cos, phi_sin | ||
|
|
||
|
|
||
|
|
@@ -121,8 +120,8 @@ | |
| provided. Further information can be found in Ruitort-Mayol et al. | ||
| drop_first: bool | ||
| Default `False`. Sometimes the first basis vector is quite "flat" and very similar to | ||
| the intercept term. When there is an intercept in the model, removing the first basis | ||
| vector may improve sampling. | ||
| parameterization: str | ||
| Whether to use `centred` or `noncentered` parameterization when multiplying the | ||
| basis by the coefficients. | ||
|
|
@@ -210,13 +209,6 @@ | |
| if parameterization not in ["centered", "noncentered"]: | ||
| raise ValueError("`parameterization` must be either 'centered' or 'noncentered'.") | ||
|
|
||
| self._drop_first = drop_first | ||
| self._m = m | ||
| self._m_star = int(np.prod(self._m)) | ||
|
|
@@ -474,11 +466,15 @@ | |
| self, | ||
| m: int, | ||
| scale: Optional[Union[float, TensorLike]] = 1.0, | ||
| drop_first=False, | ||
| *, | ||
| mean_func: Mean = Zero(), | ||
| cov_func: Periodic, | ||
| ): | ||
| arg_err_msg = ( | ||
| "`m` must be a positive integer as the `Periodic` kernel approximation is " | ||
| "only implemented for 1-dimensional case." | ||
| ) | ||
|
|
||
| if not isinstance(m, int): | ||
| raise ValueError(arg_err_msg) | ||
|
|
@@ -488,7 +484,8 @@ | |
|
|
||
| if not isinstance(cov_func, Periodic): | ||
| raise ValueError( | ||
| "`cov_func` must be an instance of a `Periodic` kernel only. Use the `scale` " | ||
| "parameter to control the variance." | ||
| ) | ||
|
|
||
| if cov_func.n_dims > 1: | ||
|
|
@@ -498,6 +495,7 @@ | |
|
|
||
| self._m = m | ||
| self.scale = scale | ||
| self.drop_first = drop_first | ||
|
|
||
| super().__init__(mean_func=mean_func, cov_func=cov_func) | ||
|
|
||
|
|
@@ -577,8 +575,7 @@ | |
| ppc = pm.sample_posterior_predictive(idata, var_names=["f"]) | ||
| """ | ||
| Xs, _ = self.cov_func._slice(Xs) | ||
| phi_cos, phi_sin = calc_basis_periodic(Xs, self.cov_func.period, self._m) | ||
| J = pt.arange(0, self._m, 1) | ||
| # rescale basis coefficients by the sqrt variance term | ||
| psd = self.scale * self.cov_func.power_spectral_density_approx(J) | ||
|
|
@@ -603,21 +600,28 @@ | |
| (phi_cos, phi_sin), psd = self.prior_linearized(X - self._X_mean) | ||
|
|
||
| m = self._m | ||
|
|
||
| if self.drop_first: | ||
|
There was a problem hiding this comment. I wonder is we would like to inform the user when There was a problem hiding this comment. Agree this would be very cool. The basis vectors are made in pytensor. I tried using |
||
| # Drop first sine term (all zeros), drop first cos term (all ones) | ||
| beta = pm.Normal(f"{name}_hsgp_coeffs_", size=(m * 2) - 2) | ||
|
There was a problem hiding this comment. Which one is the old case with There was a problem hiding this comment. Neither, now
True, need to double check this is all good, especially if someone uses an odd number There was a problem hiding this comment. I'm okay with the 
If you go against the maths, please add some comments to explain why. Otherwise pymc devs in 2034 might get confused when they look at appendix b in the original paper. There was a problem hiding this comment.
Love the optimism! There was a problem hiding this comment.
Yup it does. OK, I made a dumb mistake. When sin(x) = 0, cos(x) is either 1 or -1 (been a while since trig class), not cos(x) always equals one like I wrote. So it only happens sometimes that the cosine term is all ones. Could be all -1's or be a mix of 1's and -1's. So To see what I'm trying to explain (sorry kind of poorly), try running the code below, but change
Samples of the prior predictive for the gp look good in all cases. import pymc as pm
import numpy as np
import matplotlib.pyplot as plt
import pytensor
import pytensor.tensor as pt
X = np.linspace(0, 10, 1000)[:, None]
with pm.Model() as model:
period = 1
scale = pm.HalfNormal("scale", 10)
cov_func = pm.gp.cov.Periodic(1, period=period, ls=1.0)
m = 200
gp = pm.gp.HSGPPeriodic(m=m, scale=scale, cov_func=cov_func)
X_mean = np.mean(X, axis=0)
X = pm.MutableData("X", X)
Xs = X - X_mean
# FROM pm.gp.hsgp_approx.calc_basis_periodic(Xs, cov_func.period, m=m)
w0 = (2 * pt.pi) / period # angular frequency defining the periodicity
m1 = pt.tile(w0 * Xs, m)
m2 = pt.diag(pt.arange(0, m, 1))
mw0x = m1 @ m2
(phi_cos, phi_sin), psd = gp.prior_linearized(Xs=Xs)
plt.plot(phi_cos[0, :].eval())
plt.plot(phi_sin[0, :].eval())
beta = pm.Normal(f"beta", size=m * 2)
beta_cos = beta[:m]
beta_sin = beta[m:]
cos_term = phi_cos @ (psd * beta_cos)
sin_term = phi_sin @ (psd * beta_sin)
f = pm.Deterministic("f", cos_term + sin_term)
#prior = pm.sample_prior_predictive()
#import arviz as az
#f = az.extract(prior.prior, var_names="f")
#plt.figure(figsize=(10, 4))
#plt.plot(X.eval().flatten(), f.mean(dim="sample").data)
#plt.xticks(np.arange(10));There was a problem hiding this comment. Hmm I stared at the equation again and I can't see how So There was a problem hiding this comment.
Yuuuuuup. You're right. I plotted it wrong, so sorry for spreading my confusion around! It looks like I found that the cosine term is somtimes all ones, so the issue with the extra intercept is real. But then plotted it wrong and proceeded to draw lots of incorrect conclusions from there. What the code should do is either:
Is that right @theorashid? There was a problem hiding this comment. Sounds right to me. Maybe remove the |
||
| self._beta_cos = pt.concatenate(([0.0], beta[: m - 1])) | ||
| self._beta_sin = pt.concatenate(([0.0], beta[m - 1 :])) | ||
|
|
||
| else: | ||
| # The first eigenfunction for the sine component is zero | ||
| # and so does not contribute to the approximation. | ||
| beta = pm.Normal(f"{name}_hsgp_coeffs_", size=(m * 2) - 1) | ||
| self._beta_cos = beta[:m] | ||
| self._beta_sin = pt.concatenate(([0.0], beta[m:])) | ||
|
|
||
| cos_term = phi_cos @ (psd * self._beta_cos) | ||
| sin_term = phi_sin @ (psd * self._beta_sin) | ||
| self.f = pm.Deterministic(name, self.mean_func(X) + cos_term + sin_term, dims=dims) | ||
| return self.f | ||
|
|
||
| def _build_conditional(self, Xnew): | ||
| try: | ||
| X_mean = self._X_mean | ||
|
|
||
| except AttributeError: | ||
| raise ValueError( | ||
|
|
@@ -626,14 +630,15 @@ | |
|
|
||
| Xnew, _ = self.cov_func._slice(Xnew) | ||
|
|
||
| phi_cos, phi_sin = calc_basis_periodic(Xnew - X_mean, self.cov_func.period, self._m) | ||
| m = self._m | ||
| J = pt.arange(0, m, 1) | ||
| # rescale basis coefficients by the sqrt variance term | ||
| psd = self.scale * self.cov_func.power_spectral_density_approx(J) | ||
|
|
||
| cos_term = phi_cos @ (psd * self._beta_cos) | ||
| sin_term = phi_sin @ (psd * self._beta_sin) | ||
| return self.mean_func(Xnew) + cos_term + sin_term | ||
|
|
||
| def conditional(self, name: str, Xnew: TensorLike, dims: Optional[str] = None): # type: ignore | ||
| R""" | ||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Remember to add it to the dostrings ;) Also type:
drop_first: bool = Flase