I have no little experience with algorithms. I would like to learn how to simulate a real valued Gaussian process with correlation (e.g., exponentially correlated $\Gamma(t)=e^{-\vert t \vert/2}).$

I've tried searching for such algorithms, but I cannot determine (i) how to simulate the process over the entire line instead of over some interval (using a Fourier transform) and (ii) how to give a discrete version of the algorithm.

Any examples or links to existing examples (perhaps with a minor explanation) would be extremely helpful.


I have not found any exemple of your problem but I used the Deserno text at https://www.cmu.edu/biolphys/deserno/pdf/corr_gaussian_random.pdf to code your problem as :

 cgrn = NestList[
   E^(-1/τ) # + 
     Sqrt[1 - (E^(-1/τ))^2] RandomReal[NormalDistribution[]] &, 
      RandomReal[NormalDistribution[]], 100];

As a beginner you certainly needs some explanations : Mathematica collects nearly all the programming paradigms but it is trully efficient if you use the functional programming one. So forget the loops. Here I use Nest --- more specifically NestList which insert in a list all the calculated elements of the nested function --- i.e. : $f[f[f[...[x]]...]$. Here one has $cgrn = [f[x], f^2[x],..., f^{100}[x]]$

The best thing is not to list $cgrn$ but to plot it.


enter image description here

On must now verify if it is a gaussian process. As it is stationnary one can try to plot an histogram with more data --- 1000.

 ℋ = DistributionFitTest[cgrn, Automatic, "HypothesisTestData"];
 ℋ["TestDataTable", All]

If you know some thing about tests I let you find the result which can be verified by eyesight with

  Show[Histogram[cgrn, Automatic, "ProbabilityDensity"], 
  Plot[PDF[ℋ["FittedDistribution"], x], {x, -5, 5}, 
  PlotStyle -> Thick]]

enter image description here

Now in what concerns the autocorrelation function you have many choice to verify if it replicates what you expect.

 acf = ListCorrelate[cgrn, cgrn, {1, 1}, 0];
 ListPlot[cgrn, Filling -> Axis]

enter image description here


 ListPlot[CorrelationFunction[cgrn, {2, 50}], Filling -> Axis]

enter image description here

This last picture shows that I have made a mistake. I have tried also with the following code

 cgr[0] = RandomReal[NormalDistribution[]] ;
 cgr[n_] := cgr[n] = E^(-1/\[Tau]) cgr[n - 1] + Sqrt[1 - (E^(-1/\[Tau]))^2] RandomReal[NormalDistribution[]]
 cgrn1 := Table[cgr[i], {i, 1, 1000}];

but I encounter the same problem. Nevertheless I have decided to post my code since your question has no answer for a long time and perhaps it could help you. I hope also that an other person will see the code and will find my mistake.

I add tha of course I have not taken exactly your correlation function but even with a correction I encounter the same mistake.

| improve this answer | |
  • $\begingroup$ Thank you for your answer. How does the last picture show you made a mistake? $\endgroup$ – Iced Palmer Oct 20 '16 at 5:28
  • $\begingroup$ It is not exponentially decreasing $\endgroup$ – cyrille.piatecki Oct 20 '16 at 8:50

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.