0
$\begingroup$

Problem is related with my other question

So, generate exponentially correlated row $π‘₯_{0,0},π‘₯_{1,0}...π‘₯_{𝑛,0}$ and starting from π‘₯0,0 there will be a correlated column $π‘₯_{0,1},π‘₯_{0,2}...π‘₯_{0,n}$.

During generation we used a correlation $𝜌=𝑒π‘₯𝑝(βˆ’Ξ”πœ/πœ‰)$. And $π‘₯_{𝑗,0}=\sqrt{1βˆ’πœŒ^2}Rand[NormalDistribution[πœ‡,𝜎]]+𝜌 π‘₯_{π‘—βˆ’1,0}$

So...During 2 dimensions it should be $𝜌=𝑒π‘₯𝑝(βˆ’Ξ”πœπ‘₯/πœ‰π‘₯βˆ’Ξ”πœπ‘¦/πœ‰π‘¦)$. And in isotropic space $𝜌=𝑒π‘₯𝑝(βˆ’2Ξ”πœ/πœ‰)$.

But how to generate $π‘₯_{1,1}$ correctly with this $𝜌$? This one with square mean looks absolutely wrong: $π‘₯_{1,1}=\sqrt{1βˆ’πœŒ2}Rand[NormalDistribution[πœ‡,𝜎]]+𝜌\sqrt{π‘₯_{1,0} π‘₯_{0,1}}$

Tried to generate with arithmetic mean: $π‘₯_{1,1}=\sqrt{1βˆ’πœŒ2}Rand[NormalDistribution[πœ‡,𝜎]]+𝜌(π‘₯_{1,0}+ π‘₯_{0,1})/2$ But it look not correlated in Y direction... Any way to make it better?

my code:

ΞΌ = 0; Οƒ = 1; ΞΎ = 10; Ο„ = 1; l = {{0, 0, 
   RandomVariate[NormalDistribution[ΞΌ, Οƒ], 
     1][[1]]}}; size = 100; ρ = E^(-(
   Abs[Ο„]/ΞΎ)); ρ2 = E^(-2 Abs[Ο„]/ΞΎ);
For[j = 1, j < size, j++,
  Subscript[x, 0] = {l[[1, 3]]}; 
  Subscript[x, 
   j] = ρ2 Subscript[x, j - 1][[1]] + 
    Sqrt[1 - ρ2^2]
      RandomVariate[NormalDistribution[ΞΌ, Οƒ], 1]; 
  AppendTo[l, {j, 0, Subscript[x, j][[1]]}]];

For[i = 1, i < size, i++,
  Subscript[x, 
   0] = ρ2 l[[1 + size (i - 1), 3]] + 
    Sqrt[1 - ρ2^2]
      RandomVariate[NormalDistribution[ΞΌ, Οƒ], 1]; 
  AppendTo[l, {0, i, Subscript[x, 0][[1]]}];
  For[j = 1, j < size, j++,
   Subscript[x, 
    j] = ρ2 (Subscript[x, j - 1][[1]] + l[[size i + j, 3]])/2 + 
     Sqrt[1 - ρ2^2]
       RandomVariate[NormalDistribution[ΞΌ, Οƒ], 1]; 
   AppendTo[l, {j, i, Subscript[x, j][[1]]}]]];
ListContourPlot[l]
ListPlot3D[l]
Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ The standard way to generate correlated noise, is to: 1) Compute the autocorrelation function; 2) Fourier Transform this to get the Power spectrum; 3) Generate independent random numbers; 4) scale by the sqrt of the Power spectrum; 5) Fourier transform; 6) Verify that you have the autocorrelation function you wanted. (This will work in any number of dimensions) $\endgroup$
    – mikado
    Mar 9, 2019 at 22:24
  • 1
    $\begingroup$ @mikado: Not sure where I heard this many years ago: "We must love standards because we have so many of them." For me the standard way for this problem is to generate random samples directly from a multivariate normal with the desired covariance structure. $\endgroup$
    – JimB
    Mar 9, 2019 at 23:03
  • $\begingroup$ @JimB if I want to generate (say) a 1000x1000 image of correlated noise, my method involves generating 10^6 random samples, weighting each one with a simple function and 2D FFT (quick). Yours involves generating a 10^6 x 10^6 covariance matrix. Mathematica would struggle to generate the resulting random vector. $\endgroup$
    – mikado
    Mar 9, 2019 at 23:37
  • $\begingroup$ @mikado Yep, there would be a bit of a struggle. My weak/feeble defense is that I didn't see the OP state any size requirements. $\endgroup$
    – JimB
    Mar 9, 2019 at 23:42

1 Answer 1

Reset to default
3
$\begingroup$

As noted by @mikado generating the random sample directly is more than a bit sluggish. If you have something on the order of 10,000 points (a 100 x 100 array), then the following code is simple and straightforward but takes about 15 minutes:

n = 100;
ρ = 0.1;
x = Flatten[Table[{i, j}, {i, n}, {j, n}], 1];
Σ = Table[If[i == j, 1, Exp[-ρ Norm[x[[i]] - x[[j]]]]], {i, Length[x]}, {j, Length[x]}];
z = Transpose[{x[[All, 1]], x[[All, 2]], 
  Flatten[RandomVariate[MultinormalDistribution[ConstantArray[0, Length[x]], Ξ£], 1]]}];]
ListContourPlot[z]
ListPointPlot3D[z]

Contour plot Scatter of sample points

Improve this answer
$\endgroup$

Your Answer

By clicking β€œPost Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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