# Generate Gaussian process with squared exponential covariance function

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:

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

I believe this model also underpins the Predict function of Mathematica described here.

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

• Perhaps there are several random processes which have the same shape Covariance function. I'd be interested to know a bit more about the underlying mathematics. Dec 12, 2018 at 21:01
• Chapter 4.5 of Ripley (1987) Stochastic Simulation seems useful, though I don't fully understand it yet. Dec 13, 2018 at 12:44
• You should note when you've posted the same question simultaneously on two Stack Exchange forums: stats.stackexchange.com/questions/381729/… 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. (Although I see that one question is asking for a general approach and here you're asking for how to do it in Mathematica.)
– JimB
Dec 14, 2018 at 6:34
• @JimB True, with this question, I am specifically asking whether there is a Mathematica implementation which allows for simulation of random processes following certain requirements (such as it being stationary, Gaussian etc). Dec 15, 2018 at 15:30
• I should have noticed this before: Don't you mean $C(y_i,y_j)=\theta e^{-||x_i-x_j||^2}$? $y$ is the random variable and $x$ represents a fixed and known "location". (That's how I answered the question.)
– JimB
Jan 13, 2019 at 16:12

Edit: Previously I used σ = 1 and wrote that this approach only works for small $$n$$. But if one uses a machine precision number σ = 1., then much larger $$n$$ can be used. Using $$n=300$$ takes about 2.5 seconds to generate a sample.

If you just need a small number of distinct points in a random field, then the following brute force approach will work:

n = 300;
μ = ConstantArray[0, n];
σ = 1.;
x = Range[n];
Σ = Table[If[i == j, σ^2, σ^2 Exp[-(x[[i]] - x[[j]])^2]], {i, n}, {j, n}];
SeedRandom[12345];
y = RandomVariate[MultinormalDistribution[μ, Σ], 1]
(* {{-0.988754, 1.31061, 0.0974487,
0.0611094, -1.44568, -0.287135, -0.0662409, -0.764446, 0.197805, 0.584592,
-0.571173, 0.728825, 0.303415, -1.03507, -1.15253, -0.54701, -1.42998, -1.41337, -1.59949,
-1.87218, -0.0129242, 0.612778, 0.647016, 1.31446, 1.15284, -0.0106161,
-1.22277, 1.13532, 0.595458, -0.540409, 0.812264, -0.0404555, -1.58148,... *)


And by "small" I mean $$n \le 300$$. Also, if the correlation is positive and needs to decrease a bit faster, then you might want to consider $$\exp({-\theta (x_i-x_j)^2})$$.

• I guess $\Sigma$ gets huge quickly, as it scales with $n^2$. Though fixing distant $x$'s to zero correlation might help. No other (approximate) solutions using e.g. Fourier Transform to simulate large fields in many dimensions? Dec 13, 2018 at 18:04
• There are definitely other solutions for large $n$. I just wanted to give you a simple approach for small $n$. Others will almost certainly provide more efficient approaches. So I suggest "un-doing" the acceptance and wait at least a couple of days to accept an answer (but only if the answer really answers your question).
– JimB
Dec 13, 2018 at 18:32
• @LBogaardt Indeed, fixing small numbers to zero will help. You can use SymmetrizedArray to generate a more efficient sparse representation of the covariance matrix that does not need to store duplicate matrix entries. Dec 14, 2018 at 9:37