Add Spline Interpolation using basis functions

[ferrine]

GPs for time series or spatial data often require a lot of compute. Having dense grid to evaluate GP on, makes GP prohibitively large. Spline interpolation is a good way to reduce the GP domain maintaining coverage for the whole grid.

In the PR I aim to implement 1d spline interpolation for time series. Adding a simple util can ease the usage of GPs on large grids

Example API

import pymc as pm
import numpy as np
import pymc_experimental as pmx
import matplotlib.pyplot as plt
half_months = np.linspace(0, 365, 12*2)
with pm.Model(coords=dict(knots_time=half_months, time=np.arange(365))) as model:
    kernel = pm.gp.cov.ExpQuad(1, ls=365/12)
    # ready to define gp (a latent process over parameters)
    gp = pm.gp.gp.Latent(
        cov_func=kernel
    )
    y_knots = gp.prior("y_knots", half_months[:, None], dims="knots_time")
    y = pm.Deterministic(
        "y", 
        pmx.utils.spline.bspline_interpolation(y_knots, n=365, degree=3),
        dims="time"
    )
    trace = pm.sample_prior_predictive(1)
plt.scatter(half_months, trace.prior["y_knots"].sel(chain=0, draw=0))
plt.plot(np.arange(365), trace.prior["y"].sel(chain=0, draw=0))

[twiecki]

Did you borrow from @aseyboldt's BayesAlpha/BayesBeta?

[ferrine]

Did you borrow from @aseyboldt's BayesAlpha/BayesBeta?

I added this info into Notes

[ferrine]

I think I should also add evaluation on a non-regular grid

[ferrine]

@ricardoV94 or @twiecki do you have any other suggestions? I'm ready to just merge