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:

$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$.

|improve this question
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – LBogaardt Dec 12 '18 at 21:01
  • $\begingroup$ Chapter 4.5 of Ripley (1987) Stochastic Simulation seems useful, though I don't fully understand it yet. $\endgroup$ – LBogaardt Dec 13 '18 at 12:44
  • $\begingroup$ 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.) $\endgroup$ – JimB Dec 14 '18 at 6:34
  • $\begingroup$ @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). $\endgroup$ – LBogaardt Dec 15 '18 at 15:30
  • $\begingroup$ 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.) $\endgroup$ – JimB Jan 13 at 16:12
2
$\begingroup$

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})$.

|improve this answer
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – LBogaardt Dec 13 '18 at 18:04
  • 1
    $\begingroup$ 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). $\endgroup$ – JimB Dec 13 '18 at 18:32
  • 1
    $\begingroup$ @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. $\endgroup$ – Sjoerd Smit Dec 14 '18 at 9:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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