  We present a new strategy for learning the functional relation between a pair
of variables, while addressing inhomogeneities in the correlation structure of
the available data, by modelling the sought function as a sample function of a
non-stationary Gaussian Process (GP), that nests within itself multiple other
GPs, each of which we prove can be stationary, thereby establishing sufficiency
of two GP layers. In fact, a non-stationary kernel is envisaged, with each
hyperparameter set as dependent on the sample function drawn from the outer
non-stationary GP, such that a new sample function is drawn at every pair of
input values at which the kernel is computed. However, such a model cannot be
implemented, and we substitute this by recalling that the average effect of
drawing different sample functions from a given GP is equivalent to that of
drawing a sample function from each of a set of GPs that are rendered
different, as updated during the equilibrium stage of the undertaken inference
(via MCMC). The kernel is fully non-parametric, and it suffices to learn one
hyperparameter per layer of GP, for each dimension of the input variable. We
illustrate this new learning strategy on a real dataset.
