2
$\begingroup$

In a (stationary) Gaussian Process, values which are closeby are more similar than values far away from each other. The correlation function tends to zero as distance increases. Often, one models the decaying correlation functon $C$ as:

$C(x_i, x_j) = \theta \, e^{-||x_i - x_j||^2}$

I believe this model also underpins the Kriging method of interpolation.

However, how does one generate (i.e. simulate on a computer) a random field with such a property? You may, for simplicity, assume it's a one dimensional function $x(t)$ with mean $\mu = 0$ and standard deviation $\sigma = 1$.

$\endgroup$
1
  • $\begingroup$ You should note when you've posted the same question simultaneously on two Stack Exchange forums: mathematica.stackexchange.com/questions/187776/…. That way someone from one forum might see a response on the other forum and result in a better answer. Also, one doesn't want to waste the time of those answering if the answer is already available on one of the forums. $\endgroup$
    – JimB
    Dec 14, 2018 at 6:31

2 Answers 2

2
$\begingroup$

A Gaussian process is a probability distribution over functions, parameterized by a mean function $\mu(x)$ and covariance function $C(x, x')$. For any set of points $\{x_1, \dots, x_n\}$, the corresponding function values $y = [f(x_1), \dots, f(x_n)]^T$ have a joint Gaussian distribution with mean $m = [\mu(x_1), \dots, \mu(x_n)]^T$ and covariance matrix $K$, where $K_{ij} = C(x_i, x_j)$:

$$p(y \mid m, K) = \text{det}(2 \pi K)^{-\frac{1}{2}} \exp \left[ -\frac{1}{2} (y-m)^T K^{-1} (y-m) \right]$$

To generate a Gaussian process, you would simply pick the mean and covariance functions. To sample from this Gaussian process, you would first pick the points $\{x_1, \dots, x_n\}$ at which the function is to be evaluated. Compute $m$ and $K$ as above. Then, generate the function values by sampling vector $y$ from a Gaussian distribution with mean $m$ and covariance matrix $K$.

$\endgroup$
2
  • $\begingroup$ That makes sense. Though I can imagine this is not the most computanionally efficient method as the correlation matrix scale with the square of the number of points and sampling from it will be difficult. $\endgroup$
    – LBogaardt
    Dec 13, 2018 at 14:47
  • $\begingroup$ I was expecting it had something to do with the Fourier transform, similar to Gaussian Random Fields with a power-law spectrum. Do you know if there is any connection? $\endgroup$
    – LBogaardt
    Dec 13, 2018 at 14:48
0
$\begingroup$

user20160 is right. Just to add a little technical help: you simulate data from a Gaussian distribution with covariance matrix $K$ by calculating the Cholesky-Decomposition $K = L L^\top$ and generating a vector of independent Gaussian RVs $Z$. Then $Y = L Z$ has the wanted distribution.

In R, you can do this (or directly simulate a Gaussian process) with the help of the packages RandomFields and RandomFieldsUtils.

$\endgroup$
4
  • $\begingroup$ Great, this solves the issue of how to sample efficiently from such a complex distribution, but given $n$ sample points, matrix $K$ is still $n^2$ in size, which can be huge... wouldn't this be very inefficient in most cases? $\endgroup$
    – LBogaardt
    Dec 13, 2018 at 14:51
  • $\begingroup$ The cholesky-decomposition and the matrix multiplication $LZ$ become computational issues before the $n^2$-times evaluation of $C(h)$ should give you any trouble. However, the biggest issue for me usually is that $K$ does not fit into my RAM anymore. But it's the best solution I know $\endgroup$
    – nope
    Dec 13, 2018 at 17:29
  • $\begingroup$ I guess most elements of $K$ will be close to zero, so perhaps saving only the correlations between closeby points can reduce memory storage of $K$. It's also symmetric, and approximate rank deficient, so e.g. SVD decomposition might help. Anyway, this brute-force method still strikes me as suboptimal. $\endgroup$
    – LBogaardt
    Dec 13, 2018 at 18:08
  • $\begingroup$ If your interested in approximate solutions: If $X_i \sim N(0,1)$ i.i.d and $Y_i \sim F$ i.i.d. with $F$ being the measure needed for the representation of the covariance function $C$ with Bochner's theorem, then $Z(t) = n^{-0.5} \sum_{j=1}^n X_j \exp(i <t,Y_j>)$ is a Gaussian process with covariance function $C$ for $n \rightarrow \infty$. $\endgroup$
    – nope
    Dec 13, 2018 at 18:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.