Trying to undestand the concept of Gaussian Process Having been trying to understand the concept of Gaussian Process but I still struggle to understand the concept of Gaussian Process. Every single videos on YouTube about Gaussian Process, they all start with an axis with some indexes and another axis with function values. But what becomes the index axis and what becomes the y axis? 
If I have 3 different variables such as height, weight, and IQ. Say I want to know relationship between 2 independent variables (height and weight) and IQ. If I set IQ as real values, then I may have a linear regression model to find the relationship. Height and weight will take x-axes and IQ will take y-axis in Linear Regression but if I want to apply Gaussian Process to my Linear Regression(assuming it is appropriate to apply Gaussian Process to my Linear Regression), it looks like the only variable I care for Gaussian Process is the IQ values which is y-axis and the x-axis will be a collection of indices(I don't know what to call it). What am I not understanding and how should I think of the Y-axis and X-axis?
 A: The X axis is actually really a (n index) set. In your case (if we ignore the restrictions that there is no human being with a negative height and height is bounded and so forth) the index set is $\mathbb{R}^2$. For every possible combination of height and weight $(h,w) \in \mathbb{R}^2$ you have a new random variable $X_{(h,w)}$. The whole GP really is the collection of those infinitely many random variables GP$=(X_{(h,w)})_{(h,w) \in \mathbb{R}^2}$. The Y axis is just $\mathbb{R}$. You collect a series of heights and weights and IQs $(h_i, w_i, \text{IQ}_i)_{i=1,...,n}$.
After "training" (training looks really different from other models in GPs, that's why I put it into quotation marks) you usually seek a prediction for the IQ given a human with a specific height $\tilde{h}$ and a specific weight $\tilde{w}$. GPs are a particularly nice class of models because they natively give you a probability distribution instead of just a prediction: What you want to compute is
$$p(x_{(\tilde{h}, \tilde{w})} | x_{(h_1,w_1)}, ..., x_{(h_n, w_n)})~~~~~~~~~~~ (*)$$
or if you want a single point prediction then you could also compute the mean of that distribution:
$$\text{"mean of" } p(x_{(\tilde{h}, \tilde{w})} | x_{(h_1,w_1)}, ..., x_{(h_n, w_n)}) = \int_{\mathbb{R}} x_{(\tilde{h}, \tilde{w})} \cdot p(x_{(\tilde{h}, \tilde{w})} | x_{(h_1,w_1)}, ..., x_{(h_n, w_n)}) dx_{(\tilde{h}, \tilde{w})} = E[X_{(\tilde{h}, \tilde{w})}|X_{(h_1,w_1)}=x_{(h_1,w_1)}, ..., X_{(h_n,w_n)}=x_{(h_n,w_n)}]$$
(it is not trivial to show that these expressions really coincide)
Three further questions should maybe answered: 


*

*How does training look like? Where do we need the values of the past IQs?

*How should we ever be able to compute this weird distribution $(*)$?

*Where does the kernel come in?


-- 2. --
The assumption of a Gaussian process is that for each finite collection of indices, for example, $((\tilde{h}, \tilde{w}),(h_1,w_1), ..., (h_1,w_1))$ the random variable (or random vector, i.e. vector of random variables, to be precise) $X_{((\tilde{h}, \tilde{w}),(h_1,w_1), ..., (h_1,w_1))} := (X_{(\tilde{h}, \tilde{w})},X_{(h_1,w_1)}, ..., X_{(h_1,w_1)})$ is normally distributed with some mean and variance. That means that
$$(*) = \frac{p(x_{(\tilde{h}, \tilde{w})}, x_{(h_1,w_1)}, ..., x_{(h_n, w_n)})}{p(x_{(h_1,w_1)}, ..., x_{(h_n, w_n)})}$$
is a quotient of two normal distributions... so it is not too hopeless and actually, using a weird inversion observation from matrix theory one can derive interesting facts about this expression, for example, it is again a normal distribution and its mean and variance are connected closely to the ones from the numerator and denominator, see the formulae on wikipedia.
-- 3. --
The precise assumption is that for a general Gaussian process GP=$(X_i)_{i \in I}$ we have that for every finite subset of indices $(i_1, ..., i_n)$, the random variable $X_{i_1}, ..., X_{i_n}$ has a Gaussian distribution. In particular, we have the mean of each random variable $\mu_i = E[X_i]$ and for each pair of random variables $X_i, X_j$ we have the covariance of them $\text{Cov}(X_i, X_j)$. Since we know already that $X_{i_1}, ..., X_{i_n}$ has a Gaussian distribution and we know its parameters $\mu_{i_1}, ..., \mu_{i_n}$ and the covariance matrix, we can essentially describe all 'finite parts' of the GP with the two functions $m : I \to \mathbb{R}, m(i) = \mu_i = E[X_i]$ and $k:I \times I \to \mathbb{R}, k(i, j) = \text{Cov}(X_i, X_j)$. This is what people call mean and kernel function.
Coming back to 2. we can now compute this expression but it depends on the basic ingredients of a Gaussian process, namely the kernel function $k$ and the mean function $m$ (i.e. a complicated matrix inversion in a matrix that consists of evaluations of these functions at the different involved indices). That is why we first need to specify these and then we can actually compute the predictive distributions or the point predictions. That also covers 3.
Coming back to 1.: Where do we actually need the past IQ values? Well, in order to execute training (i.e. do matrix inversion) and then compute the predictions we need to specify $m$ and $k$. Usually people use $m=0$ and in your example, for $k$ they would probably start with something like
  $$k((w,h),(w',h'))=e^{-||(w,h)-(w',h')||^2}$$
i.e. the covariance of two points decreases negatively exponentially with their distance in $\mathbb{R}^2$. We can view that as an answer to the question of how much an IQ at a point $(w,h)$ influences the value of the IQ at a different point $(w', h')$. The next thing you probably realise is that model will not perform very well. Why? Because we do not know whether that function actually transmits the information correctly (or does it over/undershoot?). In order to resolve this we introduce a hyperparameter
$$k((w,h),(w',h'))=e^{-\gamma ||(w,h)-(w',h')||^2}$$
The IQ values in the past are actually used in order to determine the 'best' choice of this hyperparameter by setting up the likelihood function for
$$p(x_{(w_1,h_1)}, ..., x_{(w_n, h_n)})$$
(maximise this expression as a function of $\gamma$ and then use argmax as hyperparameter), where for $x_{(w_1,h_1)}, ..., x_{(w_n, h_n)}$ we use the actually observed values for these random variables, namely $x_{(w_1,h_1)}=\text{IQ}_1, ..., x_{(w_n, h_n)}=\text{IQ}_n$
