Gaussian process regression

import matplotlib.pyplot as plt
import numpy as np

from pyinla.model import *
from pyinla.spde import *
from pyinla.utils import set_pardiso_license

set_pardiso_license("/Users/js5013/.pardiso.lic")

StrVector with 1 elements.

'SUCCESS: PARDISO IS INSTALLED AND WORKING'

n_train = 250
n_test = 400
sigma_obs = 0.3

fn = (
    lambda x: 20.0 * np.sin(x / 40.0)
    + 0.1 * np.power(x, 3.0)
    + 0.1 * np.power(1.0 + x, 3.0)
    + 1.2 * np.cos(4.0 * x)
)

x = np.linspace(-1, 1, n_train)
y = fn(x)
y += sigma_obs * np.random.normal(size=(n_train,))
y -= np.mean(y)
y /= np.std(y)

x_test = np.linspace(-1.3, 1.3, n_test)

x_true = np.linspace(-1.3, 1.3, 1000)
y_true = fn(x_true)

plt.scatter(x, y)
plt.plot(x_true, y_true)

[<matplotlib.lines.Line2D at 0x38d1ff2e0>]

mesh = mesh_1d(x_test, degree=2)

sigma = 3 * np.std(y)
rho = 0.5 * np.max(x) - np.min(x)
spde = spde2_pcmatern(mesh, alpha=2, prior_range=rho, prior_sigma=sigma)

A = make_projection_matrix(mesh, loc=x)
Ap = make_projection_matrix(mesh, loc=x_test)
index = spde.make_index("ix")

st_est = inla_stack(tag="est", data=dict(y=y), A=A, effects=dict(), index=index)
st_pred = inla_stack(
    tag="pred",
    data=dict(y=np.repeat(np.nan, x_test.size)),
    A=Ap,
    effects=dict(),
    index=index,
)
st_full = combine_stacks(st_est, st_pred)

formula = "y ~ -1 + f(ix, model=spde)"
data = stack_data(st_full) | dict(spde=spde)
data["spde"] = spde
res = inla(
    formula,
    data=data,
    family="gaussian",
    control_predictor=dict(A=stack_A(st_full), compute=True),
)

res

Time used:
     = 1.83,  = 0.188,  = 0.0102,  = 2.03 
Random effects:
  Name	  Model
    ix SPDE2 model

Model hyperparameters:
                                         mean   sd 0.025quant 0.5quant
Precision for the Gaussian observations 12.42 1.16      10.26    12.38
                                        0.975quant  mode
Precision for the Gaussian observations      14.79 12.31

Marginal log-Likelihood:  -81.72 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

pred_ix = stack_index(st_full, tag="pred")
fitted = pd.DataFrame(res.get_summary("fitted.values"))
fitted = fitted.iloc[pred_ix]
fitted["x"] = x_test

plt.plot(fitted["x"], fitted["mode"])
plt.fill_between(
    fitted["x"], fitted["0.025quant"], fitted["0.975quant"], alpha=0.3, color="b"
)
plt.scatter(x, y)
plt.plot(x_true, y_true)

[<matplotlib.lines.Line2D at 0x38d305270>]