Calculating the expectation of a GPR over a given domain

[neildhir]

Hi, I am trying to calculate the expected value of a Gaussian process regression over a specific domain.

Let me start with an MWE to describe the problem.

MWE

import numpy as np
from torch import tensor
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from gpytorch.mlls import ExactMarginalLogLikelihood

# The true function of interest
f = lambda x: np.exp(-x[:, 0] ** 2 - np.sin(3.0 * x[:, 0]) ** 2) # domain=(-3, 3)

# Noisy samples of f are simulated using this
def sample_noisy_f(n, noise_std):
    x = np.linspace(-3, 3, n).reshape(-1, 1)
    noise = np.random.normal(0, noise_std, n)
    return (x, f(x) + noise)

noise_std = 0.01
x,y = sample_noisy_f(15, noise_std)    

# A GPR is fitted
def get_fitted_model(train_X:tensor, train_Y:tensor) -> SingleTaskGP:
    model = SingleTaskGP(train_X=train_X, train_Y=train_Y)
    mll = ExactMarginalLogLikelihood(likelihood=model.likelihood, model=model)
    fit_gpytorch_mll(mll) # If training succeeded then model is returned in eval() mode
    return model

model = get_fitted_model(tensor(x.reshape(-1,1)), tensor(y.reshape(-1,1)))    

Now, we are interested in:

$\mathbb{E} [ \hat{f}(x) ] = \int_{\Omega} \hat{f}(x) p(x) \text{d}x$

The relationship between $x$ and $y$ is estimated using the above GP yielding estimator $\hat{f}$ (i.e. model above). I am interested in the expectation over the domain [-3,3] so

$\mathbb{E} [ \hat{f}(x) ] = \int_{-3}^3 \hat{f}(x) p(x) \text{d}x$.

The integration itself I am not too worried about (Bayesian quadrature or MC samples I expect should be the simplest way to get hold of it). Rather, I am stuck on what I think is the much simpler item of getting hold of p(x) i.e. the weight for each output in the above regression.