# Randomness in a Gaussian Process

I'm somewhat confused about the standard setup of Gaussian Processes in the machine learning literature.

In the classic setup, we have a dataset $$D=\{\mathbf{x}, y_i \}_{i=1}^N$$ and we have that $$y=f(\mathbf{x}_i)$$ (which may or may not be corrupted by noise). I am confused exactly as to whether or not $$\mathbf{x}$$ is a random variable or not. It's usually stated that $$\mathcal{X}$$ is the index set so that, $$\{f(\mathbf{x}):\mathbf{x}\in \mathcal{X}\}$$ is a collection of random variables. From this perspective, it does not seem to make sense to think of $$\mathbf{x}$$ as a random variable. However, we also derive posterior distributions such as $$p(f|y,\mathbf{x})$$, where here we are now conditioning on $$\mathbf{x}$$, which does not make sense unless $$\mathbf{x}$$ is the realization of some random variable. How can I reconcile these two points together? Is writing the posterior distribution in this way just an abuse of notation?

The inputs $$\mathbf x \in\mathcal X$$ are not random variables, but deterministic "indices" : for instance, if $$\mathcal X =[0,T]$$, then they could represent different times, or if $$\mathcal X= \mathbb [0,1]^2$$, they could represent different positions on a map.
As you said, in Gaussian process regression, the values $$f(\mathbf x)$$ for each $$\mathbf x \in \mathcal X$$ are assumed to be random variables instead of deterministic (that is, we assume that $$\big(f(\mathbf x)\big)_{x\in\mathcal X}$$ is a Gaussian process). Therefore, the "sources of randomness" in the regression problem will be both the target $$f(\mathbf x_i)$$ and the noise $$\varepsilon(\mathbf x_i)$$ for all $$1\le i\le N$$.
With that being said, I understand your confusion about the notations commonly found in the literature. The short story is simply that when we write $$p(\mathbf y \mid\mathbf x,D)$$, the $$x$$ in the conditioning should be mostly understood as a notation to mean that we are looking at the conditional distribution of $$y$$ specifically at point $$\mathbf x$$, and not over the whole space $$\mathcal X$$.
To illustrate my point, consider the standard example : let $$\mathbf x^*\in \mathcal X$$ be a "new point" (i.e. not in the dataset), and consider the problem of finding the conditional distribution of $$f(\mathbf x^*)$$ given our observations $$D$$. Using standard notations, we denote $$f(\mathbf x^*)\equiv f^*$$ and $$\mathbf f \equiv \left[f(\mathbf x_1),\ldots,f(\mathbf x_N)\right]^T \sim \mathcal N(\boldsymbol\mu_N, \Sigma_N)$$ So now, we can write the conditional distribution of $$f^*$$ given $$D$$ by marginalizing as follows : $$P(f(\mathbf x^*)\mid D) \equiv P(f^*\mid \mathbf x^*,D) = \int P(f^*\mid \mathbf f, \mathbf x^*,D) P(\mathbf f\mid \mathbf x^*,D)\ d\mathbf f$$ Now to compute this integral, we need to express the two factors in the integrand in closed-form. For the second factor, $$\mathbf f$$ is not related to $$\mathbf x^*$$ hence it simplifies to $$P(\mathbf f\mid \mathbf x^*,D)= P(\mathbf f\mid D)$$, which may readily be found using Bayes' rule and the known distributions of $$\mathbf f$$ and $$(y_i)_{1\le i \le N}$$. For the first factor however, we have $$\begin{bmatrix}\mathbf f\\ f^*\end{bmatrix} \sim \mathcal N(\boldsymbol \mu(\mathbf x^*),\Sigma(\mathbf x^*))$$ From which $$P(f^*\mid \mathbf f, \mathbf x^*,D)$$ can be readily computed using properties of Gaussian distributions.
So, as you can see, "conditioning on $$\mathbf x^*$$" is mostly a notational device which helps keeping track of what is varying and what isn't, and can be helpful in computations, but everything would work the exact same without it. If you still have doubts, having a look at the first two pages of these notes should hopefully help clarify some of the above explanations.