Learning a Gaussian Process from function observations (not GP regression) Suppose we have a set of observations, where each observation represents a function. For example, our set is $\{f_1, f_2, ..., f_n\}$ where each $f_i = \{(x_1, y_1), (x_2, y_2), ..., (x_{p_i}, y_{p_i})\} $ representing evaluations of that function.
Now I want to learn a distribution over functions from these function observations, for example a Gaussian Process. What methods can you use for this?
Please note this is not GP regression, where you get a posterior distribution over functions representing your uncertainty over a single true function. I'm instead interested in maintaining a distribution over functions from which you sample genuinely different functions.
As a concrete example, each function observation $f_i$ could represent temperature over the course of a year, I have records from many years, and I'm interested in learning about what temperature profiles generally look like (i.e. warmer in summer, colder in winter), as well as what variations between years looks like. 
 A: First off, we need to specify a parametric form for both the mean function $\mu$ (e.g. linear function $\alpha + \beta^Tx$) and the covariance $k$ (e.g. RBF or Matérn). These parameters $\alpha$, $\beta$, $\sigma^2$ etc. are often known as hyperparameters in GP regression.
In the general case when the function paths are not i.i.d., we can index each with a variable $u$. Then the solution is to simply perform GP regression with the mean and covariance functions taking the forms $\mu(x, u)$ and $k((x, u), (x', u'))$, and then optimize the hyperparameters using e.g. maximum marginal likelihood or Bayesian inference. However, the complexity for the regression is $O(N^3)$, where N is the total number of points $(x, u)$, i.e. $N = \sum_{i=1}^n p_i$.
However, in the i.i.d. case, we can simply run a GP regression on each of the $n$ datasets $f_i$ individually, optimizing or Bayesian updating the same hyperparameters for all of these. In this case the complexity is (assuming the number of datapoints per dataset $p_i = p$ is fixed) $O(np^3)$ instead of $O(n^3p^3)$.
The resulting hyperparameter values (if using maximum likelihood) or hyperparameter distribution (if using Bayesian inference) provide the learned "distribution over functions".
Thanks @Yves for the helpful comments.
