How do we sample from a Gaussian process I have one particular question on Gaussian processes.
A Gaussian process is fully characterized by $\mu$ and $\Sigma$. However, I do not understand how can we sample a (random) function from the so defined Gaussian process.
For example, in a tutorial here, how can we generate the figures in page 6?
 A: One of the usual procedures for sampling from a multivariate Gaussian distribution is as follows.
Let $\boldsymbol{X}$ have a $n$-dimensional Gaussian distribution $N(\mu,\Sigma)$. We wish to generate a sample from $\boldsymbol{X}$.

*

*First off, you need to find a matrix $A$, such that $\Sigma = AA^T$. This is possible by something called Cholesky decomposition, and you call $A$ the Cholesky factor of $\Sigma$.


*Generate a vector $\boldsymbol{Z}=(Z_1,\ldots,Z_n)^T$ of independent, standard normal variables.


*Let $\boldsymbol{X} = \mu + A\boldsymbol{Z}$.
The $\boldsymbol{X}$ in step 3 will be the sample you are looking for.
About the Cholesky factor:
I don't know which program you use for simulating, but almost all statistical programs have finding the Cholesky factor as a build-in function, e.g. in both R and Matlab, you can use the chol() function, and SAS also has a routine for finding the $A$ matrix.
Remark: It is a prerequisite for this method to work that $\Sigma$ is positive definite (instead of just positive semi-definite), but should you happen upon a degenerate case, you can reduce the dimension of your problem to a smaller space, since the distribution corresponding to a degenerate $\Sigma$ will also be degenerate.
