I need to generate two correlated gaussian random fields. As far as I know, this question and this other provides the means to generate a single autocorrelated process. However, I am clueless about how to include a second random field, correlated with the first. The only way of doing this that I have found is on this blog and isn't a implemented in Mathematica, but in R. I am looking for a Mathematica code which does the same thing as the second blog or presents an alternative algorithm that can generate two processes as the answer to the first or second links but correlated as per the third?
$\begingroup$
$\endgroup$
1
-
$\begingroup$ What exactly do you mean by the two fields being correlated? For example, if you generate one random field, you can then generate a second field just by adding noise to the first field. Or is the correlation structure you seek more complicated than that? $\endgroup$– JimBAug 5, 2021 at 19:56
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
2
I'm not sure what you mean by "correlated" in "correlated gaussian random fields" but you can generate two independent fields and then compare the first with the sum of the first and second. The individual points are then correlated.
Stealing from the excellent answer you posted, one can construct two independent fields:
GaussianRandomField[size : (_Integer?Positive) : 256,
dim : (_Integer?Positive) : 2, Pk_ : Function[k, k^-3]] :=
Module[{Pkn, fftIndgen, noise, amplitude, s2},
Pkn = Compile[{{vec, _Real, 1}},
With[{nrm = Norm[vec]}, If[nrm == 0, 0, Sqrt[Pk[nrm]]]],
CompilationOptions -> {"InlineExternalDefinitions" -> True}];
s2 = Quotient[size, 2];
fftIndgen = ArrayPad[Range[0, s2], {0, s2 - 1}, "ReflectedNegation"];
noise = Fourier[RandomVariate[NormalDistribution[], ConstantArray[size, dim]]];
amplitude = Outer[Pkn[{##}] &, Sequence @@ ConstantArray[N@fftIndgen, dim]];
InverseFourier[noise*amplitude]]
BlockRandom[SeedRandom[42, Method -> "Legacy"]; (*for reproducibility*)
grf1 = Chop[GaussianRandomField[]]];
BlockRandom[SeedRandom[12345, Method -> "Legacy"]; (*for reproducibility*)
grf2 = Chop[GaussianRandomField[]]];
GraphicsRow[{MatrixPlot[grf1, PlotLabel -> Style["1", Bold, 18]],
MatrixPlot[grf2, PlotLabel -> Style["2", Bold, 18]]}]
Now sum the two fields:
MatrixPlot[grf1 + grf2, PlotLabel -> Style["1+2", Bold, 18]]
The individual points of 1 and 1+2 are positively correlated:
points = Transpose[{Flatten[grf1], Flatten[grf1 + grf2]}];
Correlation[points][[1, 2]]
(* 0.712464 *)
ListPlot[points, PlotLabel -> Style["1 vs 1+2", Bold, 18]]
-
$\begingroup$ Thanks for the idea. I will study it. I am working with this other article in which a covariance kernel for two processes is generated and then, the realizations of said process generated using Karhunen-Loève can be split in two. I need to do that but using 2D fields instead of 1D gaussian processes. I hope I have explained myself. $\endgroup$ Aug 6, 2021 at 7:11
-
$\begingroup$ Thanks for the accept but if someone else posts something that's more relevant to your particular objective, don't hesitate to accept that answer. $\endgroup$– JimBAug 6, 2021 at 20:00