Why do sampled functions from a Gaussian Process perfectly interpolate the data? I've been studying GPs with Rasmussens book as well as a few of my old probability favorites.
I am confused by the results we're getting (doing also a lot of practice) when sampling from our Gaussian Process.
Conditional on the data $ D_n = [{X_n, Y_n}] $ for $n$ datapoints the predictive posterior distribution for a new datapoint $x$  is $GP(x) = Y(x) \sim N(\mu(x), \Sigma(x))$ Where $\mu(x)$ is linearly dependent on the $y_n$.
I do find most of it intuitive, except for the "knots" see the screenshot from Rasmussen's book
.
At first, I thought that it was just illustrative but I was able to get the same "shapes" (though not as pretty as my plotting skills are limited) with arbitrary data, though also not consistently, just most of the time, at times I don't get the knots everywhere, especially with many points.
If I understand correctly since the covariance at those points is 1, for any realization, at those points $Y(x_1) = Y(x_2)$ if $x_1 = x_2$ But only for that sampled function.
Why are all the sampled functions getting the same value?
I do understand that this intuitively is the point to the whole exercise.
We sample from possible functions in the $GP$ that map $x$, to $y$. I am just confused as on how to see it from the math.
This is probably quite the silly question and I am missing something here, I apologize in advance and thank you for your patience!

Short answer: The reason that the sampled functions from the predictive posterior distribution perfectly interpolate the data is that they are sampled from a distribution that is conditioned on the data points. In other words, the data points are not just data points, they are also points that the function must go through.
Let's think about what we're doing here. We have a bunch of data points (x_i, y_i) and we want to use a Gaussian process to interpolate between them.
First, we need to choose a kernel function. This is a function that gives the covariance between any two points. For example, the RBF kernel is given by
k(x, x') = exp(-||x - x'||^2 / (2 * sigma^2))
where ||x - x'|| is the Euclidean distance between x and x'. sigma is a parameter that controls the "smoothness" of the function.
Next, we need to choose a mean function. This is a function that gives the expected value of the function at any point. For example, we could choose the mean function to be the constant function m(x) = 0.
Now, we can use these to define a Gaussian process. This is a collection of random variables, any one of which is a function from x to y, that have the following properties:

*

*The mean function of the Gaussian process is m(x).

*The covariance function of the Gaussian process is k(x, x').

So, what does this all mean? It means that, for any two points x and x', the expected value of the function at those points is m(x) and m(x'), and the covariance between the function values at those points is k(x, x').
Now, let's think about what happens when we sample from a Gaussian process. We're essentially choosing a random function from the set of all possible functions that could have generated our data. So, if we have two points x and x' that are close together, then the function we sample is likely to have values that are close together at those points.
This is why the sampled functions from a Gaussian process perfectly interpolate the data. Because the data points are close together, the sampled functions are likely to have values that are close together at those points. And because the mean function is the same at all points, the sampled functions will have the same value at those points.
I hope that makes sense!
