Meaning of "Distribution over functions" in Gaussian Processes? I'm reviewing PyMC3's Gaussian Process documents and it's illuminated that I might have a flawed understanding of what "distribution over functions" actually means. Consider the below code:
# A one dimensional column vector of inputs.
X = np.linspace(0, 1, 10)[:, None]

with pm.Model() as model:
    # Specify the covariance function.
    cov_func = pm.gp.cov.ExpQuad(1, ls=0.1)

    # Specify the GP.  The default mean function is `Zero`.
    gp = pm.gp.Marginal(cov_func=cov_func)

    # Place a GP prior over the function f.
    sigma = pm.HalfCauchy("sigma", beta=3)
    y_ = gp.marginal_likelihood("y", X=X, y=y, noise=sigma)

...

# After fitting or sampling, specify the distribution
# at new points with .conditional
Xnew = np.linspace(-1, 2, 50)[:, None]

with model:
    fcond = gp.conditional("fcond", Xnew=Xnew)

As you can see, the (RBF) kernel hyperparameters are fixed; no priors were attached to either in pm.gp.cov.ExpQuad(1, ls=0.1). It was my understanding that Kernel(sigma_f=1, length=0.1) is a single function, not a distribution over functions.
But perhaps I'm thinking about it wrong; rather a Kernel with fixed hyperparameters generates a covariance function (often the mean function is assumed to be all zeros.) When the mean and covariance functions are supplied as parameters to multivariate gaussian, which can be sampled from, each draw is a sample from the distribution over functions.
Any thoughts?
 A: I will assume in this answer that the notion of a multivariate distribution makes sense to you. Given this, let us first consider what a "distribution of functions" means if the functions have finite domain. For example, consider the set of functions
$$f: \{1,2,\ldots,n\}\to\mathbb R$$
To specify a function $f$ in this set, we just need to define $f(1), f(2),\ldots, f(n)$. But this is equivalent to just specifying a vector in $\mathbb R^n$. To specify a distribution over functions $f: \{1,2,\ldots,n\}\to\mathbb R$, we may therefore equivalently think about this distribution over vectors in $\mathbb R^n$, which, to reiterate, I am assuming is a notion you are already familiar with. If we make the further restriction that our distribution of functions is "Gaussian" and mean 0 in the sense that $(f(1),\ldots, f(n))$ is multivariate normally distributed and mean 0, then it is a well known fact that specifying the variance-covariance matrix of this distribution suffices to specify the distribution.
Now consider attempting to define a distribution over functions $f: [a,b]\to \mathbb R$. Because $[a,b]$ is uncountably infinite, it is difficult to proceed exactly as above and define the entire distribution all at once. But what we can do instead is to try to define a way to describe any finite dimensional marginal distribution. More precisely, letting, $a \leq x_1 < x_2 < \cdots < x_n \leq b$ for arbitrary $n$ and arbitrary $x_i$, we call the finite dimensional marginal distribution at $x_1,\ldots,x_n$ to be the distribution of the vector $(f(x_1),\ldots, f(x_n))$ (where to be clear, the randomness is coming from the randomness in $f$). It is a well known result in stochastic process theory that under a condition known as tightness, the finite dimensional marginal distributions are sufficient to characterize the distribution of $f$ (i.e. if I give you a way to describe any finite dimensional marginal distribution, then there is exactly one distribution on $f$ consistent with these finite dimensional marginal distributions).
Let us finally relate this back to Gaussian processes and why specifying a kernel suffices for specifying a distribution over functions. By definition, a Gaussian process is just a distribution over functions such that all finite dimensional marginal distributions are Gaussian. If we make the assumption (as is made by default in the Gaussian process code you cite) that our distribution is mean 0 in the sense that $E[f(x)] = 0$ for all $x$, then to specify the finite dimensional marginal distributions, it suffices to simply specify the covariance between $f(x_1), f(x_2)$ for any $x_1, x_2$. But by definition, the meaning of the kernel function is precisely that $k(x_1,x_2) = \mathbb C\mathrm{ov}(f(x_1), f(x_2))$. Then, using the fact that the finite dimensional marginal distributions suffice to characterize the distribution of functions, specifying this kernel suffices to unambiguously describe a certain distribution of functions. Thus, the kernel function is a representation of a distribution over functions (in particular, a Gaussian process) in much the same way that the variance-covariance matrix is a representation of some mean 0 multivariate normal distribution.
A: I can't figure out what PyMC3 does and I find their documentation for Gaussian processes very poorly written.
But the idea of Gaussian process (GP) regression is very simple. We have
$$
y_i = f(x_i) + \sigma w_i, \quad i =1,\dots,n,
$$
and we put a GP prior on $f$, which is a distribution over function $f \sim GP(0, K(\cdot, \cdot))$. I am assuming a zero-mean process for simplicty. Here $K(\cdot,\cdot)$ is the kernel function, that is, a symmetric (positive semidefinite) bivariate function specifying the covariance matrix of a $(f(x_1),\dots,f(x_n))$ for any given collection of points $\{x_i\}$.
To simplify notation, I will write $f(x) \in \mathbb R^n$ for a vector $x = (x_1,\dots,x_n)$ to mean $f(x):= (f(x_1),f(x_2),\dots,f(x_n))$. What the above means is that
$$
\text{cov}(f(x)) = K(x,x).
$$
where $K(x,z)$ is my notation for the matrix with entries $K(x_i, z_j)$ for any vectors $x$ and $z$.
The GP prior implies that for any two vectors $x$ and $x^*$, the vectors $f(x)$ and $f(x^*)$ are jointly Gaussian as follows
$$
\begin{bmatrix}
f(x) \\ f(x^*) 
\end{bmatrix} \sim N \Big( 
\begin{bmatrix}
0 \\ 0
\end{bmatrix},
\begin{bmatrix}
K(x,x) & K(x,x^*)  \\
K(x^*,x) & K(x^*, x^*)
\end{bmatrix}
\Big).
$$
The covariance function $K(\cdot, \cdot)$ allows you to compute the joint distribution of $(f(x),f(x^*))$ for any collection of points $x$ and $x^*$.
Our original model in vector form can be written as
$
y = f(x) + \sigma w
$
and we have
$$
\begin{bmatrix}
y \\ f(x^*) 
\end{bmatrix} \sim N \Big( 
\begin{bmatrix}
0 \\ 0
\end{bmatrix},
\begin{bmatrix}
\sigma^2 I + K(x,x) & K(x,x^*)  \\
K(x^*,x) & K(x^*, x^*)
\end{bmatrix}
\Big)
$$
assuming $w$ and $f$ are independent. Assuming $\sigma^2 = 0$ for simplicity, by the standard results on the conditional distributions for Gaussians, we have that $f(x^*)$ given $y$ is distributed as a Gaussian:
\begin{align*}
f(x^*) \mid y \sim N\Big(K(x^*,x)K(x,x)^{-1} y, \;K(x^*,x^*) - K(x^*,x) K(x,x)^{-1} K(x,x^*)\Big).
\end{align*}
This is the posterior distribution of $f(x^*)$ given $y$ for a given set of points $x^*$. There is also a posterior distribution on the entire $f$ given $y$ which in practice only accessed via its finite-dimensional sections, i.e., vectors like $f(x^*)$.

In PyMC3, the line
import pymc3 as pm
cov_func = pm.gp.cov.ExpQuad(1, ls=0.1)

defines the kernel function $K(\cdot, \cdot)$. You can apply it to sample form the prior GP + noise as follows
import numpy as np
import matplotlib.pyplot as plt

n = 100
x = np.linspace(-1, 1, n)[:, None]

cov_func = pm.gp.cov.ExpQuad(1, ls=0.1)
K = cov_func(x).eval()

f = np.random.multivariate_normal(np.zeros(K.shape[0]), K, 1).T
y = f + np.random.normal(0, 0.5, (n,1))

plt.figure(figsize=(14, 4))
plt.plot(x, f,'b')
plt.plot(x, y,'ro')


That is, the code samples from a $f \sim GP(0, K(\cdot,\cdot))$ evaluated at points x to get vector f $= f(x)$ above. In other words, it samples $f(x) \sim N(0, K(x,x))$ where $K(x,x)$ is the kernel matrix corresponding to points $x = (x_1,\dots,x_n)$ (i.e., K in the above code). We are also sampling the noising observations $y = f(x) + \sigma w$ where $\sigma^2 = 0.5$.
You should now be able to use the above expressions to sample from the posterior $f(x^*) \mid y$ using your choice of new points $x^*$.
