Reproduce example from NIPS 2017 paper

[j-faria]

Hi,
I was wondering if there is an example or tutorial about using mogptk to reproduce the synthetic example shown in Fig 2 of the Spectral Mixture Kernels for Multi-Output Gaussian Processes paper. This is the figure

I'm almost sure all the relevant methods are already implemented in the MOSM model, just wondering if there is a script that could be adapted more easily to reproduce this analysis. In particular, I'm still looking for an easy way to plot the auto and cross covariances of the model.

Thanks!

[tdewolff]

Hi João, that paper was published before we created the toolkit. @felipe-tobar or @mrlz do you know more?

[j-faria]

@tdewolff, thank you for your answer, I'm getting back to this now.

Do you have any suggestions on how to plot the auto and cross covariances of the MOSM models.
Is this what the plot_cross_spectrum is showing?

mogptk/mogptk/models/mosm.py

Line 148 in 3927103

def plot_cross_spectrum(self, title=None, figsize=(12,12)):

I'm ultimately trying to estimate (and visualize) the delay between two signals (the value of ~2 in the plot above). I still didn't understand how to do this within the MOSM models.

[tdewolff]

Sorry for the late response. The cross spectrum plot is showing the PSDs of each channel and the covariances between the channels. What you are looking for is probably the kernel evaluated for a given tau between two channels (or the auto-covariance). There is no easy helper function for that, so you'd have to get your hands dirty with the core GPR library. Something along these lines works (we'll be implementing this in the library soon):

dataset = ...
model = mogptk.MOSM(dataset, Q=3)

n = 101
channel = np.ones((n,1))
fig, ax = plt.subplots(len(dataset), len(dataset), squeeze=False)
for j in range(len(dataset)):
    for i in range(len(dataset)):
        if j < i:
            continue

        x_range = dataset[j].X[0].transformed.max() - dataset[j].X[0].transformed.min()
        x_range /= 4
        tau = np.linspace(-x_range, x_range, num=n).reshape(-1,1)

        X0 = np.concatenate((j*channel,tau), axis=1)
        X1 = np.array([[i,0.0]])
        k = model.kernel(X0,X1).detach().numpy()
        ax[j,i].plot(tau, k)
plt.show()

EDIT: I've added support for this in the library with the function model.plot_kernel(), please let me know if that works for you!

[tdewolff]

Please use the model.plot_kernel() for this, closing issue