In https://stats.stackexchange.com/a/103533/83252 @juampa writes
[...] in models like Naive Bayes or Gaussian Processes you would integrate out missing variables, and choose the best option with the remaining variables.
What does (s)he mean by that? How do I use Gaussian Processing Regression on incomplete datasets, "integrating out" the missing variables?
I am mostly interested in a detailed (hopefully step by step, easy to understand) mathematical explanation.
The linked question is discussing data imputation for the purposes of building a predictive model. What I believe the accepted answer is referring to is using a Gaussian Process as a model for the missing data, conditional on the observed data.
"Integrating out" these missing variables then means marginalising the predictions of the resulting predictive model (the model built using these "imputed" data points) over the distribution of possible values the missing points could take.
I am not confident enough to give you a rigorous mathematical explanation, but I will attempt a verbose, intuitive example since this has been open for a while, and perhaps someone else can expand on it or correct me if I also misunderstand.
Suppose we wish to build a predictive model of some function $f(x,y,z)$, but the process by which we collect data on variable $z$ is messy. So sometimes we will only have access to $x$ and $y$.
In order to build our predictive model, we require $z$. One approach is to attempt to "impute" $z$, given it's previous observations and any covariance structure we believe might be present between $x$, $y$ and $z$: this is where we can use a Gaussian Process.
In this case, we suppose a Gaussian Process models $z$ as a function of $x$ and $y$. To make things a bit easier to read, assume we concatenate $x$ and $y$ into a vector, $r$ (so $r:=[x, y]$). We consequently end up with a set of different $r$ corresponding to all the data points we have collected.
For some of these datapoints, there will not be a corresponding measurement $z$ (as it was corrupted, is missing, or otherwise unavailable). If we refer to these points ($x$ and $y$ measurements missing a particular $z$) as $r^*$, and the missing values as $z^*$, our task is to compute $z^*$ at $r^*$ conditional on $r$ and $z$.
Using a Gaussian Process we can do this as follows:
Let $k$ be some appropriate covariance kernel and $\theta$ its hyperparameters, and let $\mathbf{z}$ be the stacked vector of observed $z$ (so those $z$ measurements we do have). Additionally, I will use the somewhat abusive notation of $k(a,b)$ to represent the covariance matrix obtained by evaluating $k$ pairwise between the elements of $a$ and $b$.
Note: If you're unfamiliar with how to pick a kernel or how to identify the hyperparameters, I suggest reading some of the introductory material on Gaussian Processes - I won't cover it here in the interest of brevity. I personally like this three part (1,2,3) series by Michael Betancourt.
The predictive distribution of $z^*$ is a multivariate normal distribution conditional on the currently observed values $z$, the points at which they are observed, $r$, the kernel $k$, and the kernel hyperparameters $\theta$:
$p(z^* | r, z, k, \theta) \sim \mathrm{MultiNormal} (\mu_p, \Sigma_p)$
where the "predictive" (or conditional) mean distribution $\mu_p$ is given by:
$\mu_p = k(r^*,r) [k(r, r)]^{-1} \mathbf{z}$
And the predictive covariance matrix $\Sigma$ is:
$\Sigma = k(r^*, r^*) - k(r^*, r) [k(r, r)]^{-1} k(r,r^*)$
Given this model for $p(z^* | r, z, k, \theta)$, we are able to assess the probability of $z^*$ given the observed data and the assumptions we used to select $k$.
Now, suppose we want to use these imputed values of $z^*$ in a predictive model for some other quantity $f(x,y,z)$: It is here we would "integrate out" missing variables.
Since our predictions of $z^*$ are a distribution (a normal distribution, $p(z^*)$, given we have used a Gaussian Process as a model for $z$), we probably want to capture uncertainty in the fact we do not know $z$ at these locations exactly.
Very simply, if we want to make a prediction of $f(x,y,z)$ using $x, y$ and $p(z^*)$, we should account for all the values $z^*$ could be when making this prediction. This amounts to approaching the following integral:
$p(f(x,y,z^*)) = \displaystyle\int f(x,y,z^*) p(z^*) dz^*$
This is usually approached numerically since $f$ rarely permits an analytical treatment (but it might). Typical approaches are quadrature or Monte-Carlo simulation.
The important point is that predictions made using $z^*$ are done so over the distribution of possible $z^*$. This implies these predictions are themselves a distribution, and should be assessed as such to evaluate predictive performance in a robust way,
