docs: write documentation for robustgasp module

docs/source/emulator/index.rst

@@ -20,7 +20,7 @@ Currently implemented classes in this module are:
   multiple values as output. 
 
 .. toctree::
-   :maxdepth: 1
+   :maxdepth: 2
    :hidden:
    :caption: Emulator
 

docs/source/emulator/robustgasp.rst

@@ -2,13 +2,18 @@ RobustGaSP
 ==========
 
 `RobustGaSP`, standing for `Robust Gaussian Stochastic Process Emulation`, is a
-`R` package developed by :cite:t:`Gu2019`. In this module,
-we have implemented class :class:`.ScalarGaSP` and  class :class:`.PPGaSP` 
-based on `RobustGaSP` and `rpy2`, which allows us to use `RobustGaSP` directly
-from within `Python`.
+`R` package developed by :cite:t:`Gu2019`. It implements robust methods to
+estimate the unknown parameters which can greatly improve predictive performance
+of the built emulator. In addition, it implements the `PP GaSP` (parallel partial
+Gaussian stochastic process) emulator for simulators which return multiple values
+as output, see :cite:t:`Gu2016` for more detail.
 
-Theory of Gaussian process emulation
-------------------------------------
+In :py:mod:`PSimPy.emulator.robustgasp`, we have implemented class
+:class:`.ScalarGaSP` and  class :class:`.PPGaSP` based on `RobustGaSP` and
+`rpy2`, which allows us to use `RobustGaSP` directly from within `Python`.
+
+Theory recap of GP emulation
+----------------------------
 
 A simulator represents a mapping from an input space to an output space.
 Gaussian process emulation treats a computationally expensive simulator as an
@@ -19,8 +24,8 @@ can then be used for computationally demanding analyses :cite:p:`Zhao2021a`. A r
 of Gaussian process emulation is provided below, which is adapted from `Chapter 3`
 of :cite:t:`Zhao2021b`.
 
-Gaussian process prior
-^^^^^^^^^^^^^^^^^^^^^^
+GP prior
+^^^^^^^^
 
 Let :math:`f(\cdot)` denote a simulator that one wants to approximate. It
 defines a mapping from a :math:`p`-dimensional input
@@ -71,8 +76,8 @@ dimension. The unknowns in the covariance function
 :math:`\sigma^2` and the :math:`p` correlation length parameters
 :math:`\boldsymbol{\psi}=(\psi_1,\ldots,\psi_p)^T`.
 
-Gaussian process posterior
-^^^^^^^^^^^^^^^^^^^^^^^^^^
+GP posterior
+^^^^^^^^^^^^
 
 Equations :eq:`GP` to :eq:`GPcorrelation` encode prior knowledge of the simulator
 :math:`f(\cdot)`. From the Bayesian viewpoint, this prior knowledge can be updated
@@ -133,13 +138,44 @@ updated Gaussian process, given the unknown regression parameters
 parameters :math:`\boldsymbol{\psi}`. 
 
 For any new input :math:`\mathbf{x}^*`, the output :math:`y^*` can be almost
-instantaneously predicted by the mean :math:`m'`, namely equation :eq:`conditionalGaussianM`.
-Moreover, the uncertainty associated with the prediction can be estimated by
-for example the variance :math:`\sigma^2c'`. Depending on how the unknown parameters
-(:math:`\boldsymbol{\beta}`, :math:`\sigma^2`, and :math:`\boldsymbol{\psi}`)
+instantaneously predicted, given by the mean :math:`m'`, namely equation
+:eq:`conditionalGaussianM`. Moreover, the uncertainty associated with the prediction
+can be estimated by for example the variance :math:`\sigma^2c'`. Depending on how
+the unknown parameters (:math:`\boldsymbol{\beta}`, :math:`\sigma^2`, and :math:`\boldsymbol{\psi}`)
 are learned based on the training data, the posterior can have various forms, see
 `Section 3.3` of :cite:t:`Zhao2021b`.
 
+Sample from GP
+^^^^^^^^^^^^^^
+A Gaussian process, either the prior or the posterior, defines a distribution over
+functions. Once its mean function :math:`m(\cdot)` and covariance function
+:math:`C(\cdot,\cdot)` are specified, the Gaussian process can be evaluated at a
+finite number :math:`s` of input points :math:`\{\mathbf{x}_i\}_{i=1}^{s}`.
+The joint distribution of :math:`\{y_i=f(\mathbf{x}_i)\}_{i=1}^{s}` follows a
+:math:`s`-variate Gaussian distribution. A sample of the :math:`s`-variate Gaussian
+distribution, denoted as :math:`\tilde{\mathbf{y}}:=(\tilde{y}_1,\ldots,\tilde{y}_s)^T`,
+can be viewed as a sample of the Gaussian process (evaluated at a finite number of points).
+To generate :math:`\tilde{\mathbf{y}}`, the :math:`s \times s` covariance matrix
+:math:`\boldsymbol{\Sigma}` needs to be decomposed into the product of a lower triangular
+matrix :math:`\mathbf{L}` and its conjugate transpose :math:`\mathbf{L}^T` using the Cholesky
+decomposition
+
+.. math:: \boldsymbol{\Sigma}=\left(\begin{array}{cccc}
+    C\left(\mathbf{x}_{1}, \mathbf{x}_{1}\right) & C\left(\mathbf{x}_{1}, \mathbf{x}_{2}\right) & \ldots & C\left(\mathbf{x}_{1}, \mathbf{x}_{s}\right) \\
+    C\left(\mathbf{x}_{2}, \mathbf{x}_{1}\right) & C\left(\mathbf{x}_{2}, \mathbf{x}_{2}\right) & \ldots & C\left(\mathbf{x}_{2}, \mathbf{x}_{s}\right) \\
+    \vdots & \vdots & \ddots & \vdots \\
+    C\left(\mathbf{x}_{s}, \mathbf{x}_{1}\right) & \ldots & \ldots & C\left(\mathbf{x}_{s}, \mathbf{x}_{s}\right)
+    \end{array}\right)=\mathbf{L}\mathbf{L}^T.
+    :label: Choleskydecomposition
+
+Then a sample :math:`\tilde{\mathbf{y}}` can be obtained as follows:
+
+.. math:: \tilde{\mathbf{y}} = (m(\mathbf{x}_i),\ldots,m(\mathbf{x}_s))^T + \mathbf{L}\mathbf{w},
+    :label: GPsample
+
+where :math:`\mathbf{w}:=(w_1,\dots,w_s)^T` denotes a :math:`s`-dimensional random vector
+consisting of :math:`s` independent standard normal random variables :math:`w_i,i \in \{1,\ldots,s\}`.
+
 
 ScalarGaSP Class
 ----------------