Heteroskedastic Marginal Likelihood GP

Questions
1

Is it possible to construct a heteroskedastic GP with the “Marginal Likelihood” implementation? There’s a note in the example stating “a covariance function for the noise can be given” and providing an ExpQuad+WhiteNoise kernel as an example. This does run, but the result is not heteroskedastic colab notebook.

I think I see why this doesn’t work as-is, but if its goal isn’t to achieve heteroskedasticity I’m not sure what the point of an ExpQuad+WhiteNoise noise kernel would be. Maybe another example would be more clear?

2 Likes
2

Yes it definitely is possible. I see what you’re saying with the example though, it should definitely be clearer. I think the goal there is just to demo the ExpQuad + Noise as a stand-in for any arbitrary covariance function representing noise, to show that you don’t need to just use a scalar parameter sigma. You could define your own heteroskedastic covariance and give it to the noise argument.

1 Like
3

Here’s maybe a better example with a custom noise kernel:

import numpy as np
import pymc3 as pm
import theano.tensor as tt

class Heteroskedastic(pm.gp.cov.Covariance):
    ## adapted from pm.gp.cov.WhiteNoise
    def __init__(self, sigmas):
        super().__init__(1, None)
        self.sigmas = tt.as_tensor_variable(sigmas)
        
    def diag(self, X):
        return self.sigmas
    
    def full(self, X, Xs=None):
        if Xs is None:
            return tt.diag(self.diag(X))
        else:
            return tt.alloc(0.0, X.shape[0], Xs.shape[0])

## underlying process
x = np.linspace(0, 10, 100)
f = np.sin(2 * np.pi * x * 0.1)

## heteroskedastic noise
sigma = 0.15 * x + 0.25
noise = sigma * np.random.randn(100)
y = f + noise

plt.plot(x, f)
plt.plot(x, y, '.k');

image
image722×480 44.5 KB 

## fit model
with pm.Model() as model:
    eta = pm.HalfNormal('eta', sd=5)
    ell = 3
    cov = eta**2 * pm.gp.cov.ExpQuad(1, ls=ell)
    gp = pm.gp.Marginal(cov_func=cov)
    
    a = pm.Normal('a', mu=0.0, sd=5.0)
    b = pm.HalfNormal('b', sd=5.0)
    sigmas = pm.Deterministic('sigmas', a*x + b)
    noise = Heteroskedastic(sigmas)
    gp.marginal_likelihood("y", X=x[:, None], y=y, noise=noise)
    
    tr = pm.sample(cores=1, chains=1)

## show fit
with model:
    f = gp.conditional('f', Xnew=x[:, None])
    ppc = pm.sample_posterior_predictive(trace=tr, var_names=['f', 'y'])

ub = np.percentile(ppc['f'], 95, axis=0)
lb = np.percentile(ppc['f'], 5, axis=0)

plt.plot(x, y, '.k')
plt.fill_between(x, lb, ub, color='slateblue', alpha=0.5)

image

4 Likes
4

I see, so you can do parametric heteroskedasticity with Marginal, but you need Latent for non-parametric heteroskedasticity.

Would you mind if I incorporated this approach into my example notebook, and potentially into Gumbi?

1 Like
5

Ah okay, @BioGoertz that’s exactly right. No please do, glad that it’s helpful!

1 Like