.. _grf:

Efficient sampling from Gaussian Random Fields (GRFs)
======================================================================

A GRF is a random function defined by its power spectral density (PSD)
$\hat{C}(k)$ as a function of wavevector $k$.
It is thus stationary, ie, its statistical properties are translationally
invariant.
It is also known as a Gaussian process prior.
Given a GRF model,
it is important to be able to 1) draw independent samples from the GRF
at a given set of target points, 
and 2) compute the log likelihood (under the model) of a given function values at a given set of target points.
We focus on the first application. (For the second, see reference at the
bottom of this page.)

For certain PSDs, independent samples may be generated by solving
a PDE with iid random forcing; however this introduces boundary effects.
Here we simply use Fourier methods to sample directly from any PSD
that is a smooth function of $k$,
and evaluate at arbitrary target points.
The idea is that a *regular* Fourier grid $k_j$ be used, with
spacing sufficiently small based upon the size of the spatial domain
in which target points lie. This serves as a quadrature scheme.
The GRF model corresponds to random iid coefficients
at these Fourier values, withzero mean and variance $\hat{C}(k_j)$.
These coefficients are then evaluated at target points using the
type 2 NUFFT.

The `code <https://github.com/flatironinstitute/finufft/blob/master/tutorial/samplegrf1d.m>`_ is rudimentary but commented to be self-explanatory:

.. literalinclude:: ../../tutorial/samplegrf1d.m
   :language: matlab

The result is as follows, showing 5 samples from a Yukawa or
Matern kernel, plus the pointwise variance estimated from 300 samples:

.. image:: ../pics/grf1d.png
   :width: 70%

The calculation takes around 1 second.

If $\hat{C}(k)$ has singularities or discontinuities, a different
$k$-space quadrature scheme would be needed.
For an example of this, sampling from a 2D GRF, see
`random plane waves <https://users.flatironinstitute.org/~ahb/rpws/>`_.

Further reading
~~~~~~~~~~~~~~~~

For background on Gaussian random fields, aka, Gaussian processes,
see:

  * C E Rasmussen & C K I Williams, *Gaussian Processes for Machine Learning*, the MIT Press, 2006. http://www.GaussianProcess.org/gpml
    
  * Fast Direct Methods for Gaussian Processes,
    Sivaram Ambikasaran, Daniel Foreman-Mackey, Leslie Greengard,
    David W. Hogg, and Michael O’Neil. arxiv:1403.6015 (2015)