diff --git a/src/functions-reference/deprecated_functions.Rmd b/src/functions-reference/deprecated_functions.Rmd
index 248ef599b..f51dc1032 100644
--- a/src/functions-reference/deprecated_functions.Rmd
+++ b/src/functions-reference/deprecated_functions.Rmd
@@ -9,7 +9,13 @@ if (knitr::is_html_output()) {
}
```
-## integrate_ode_rk45, integrate_ode_adams, integrate_ode_bdf ODE integrators {#functions-old-ode-solver}
+```{r results='asis', echo=FALSE}
+if (knitr::is_html_output()) {
+ cat(' * cov_exp_quad\n')
+}
+```
+
+## integrate_ode_rk45, integrate_ode_adams, integrate_ode_bdf ODE Integrators {#functions-old-ode-solver}
These ODE integrator functions have been replaced by those described in:
@@ -24,7 +30,7 @@ if (knitr::is_html_output()) {
A system of ODEs is specified as an ordinary function in Stan within
the functions block. The ODE system function must have this function
signature:
-
+
```stan
array[] real ode(real time, array[] real state, array[] real theta,
array[] real x_r, array[] int x_i);
@@ -117,7 +123,7 @@ are as follows.
* *`ode`*: function literal referring to a function specifying the
system of differential equations with signature:
-
+
```
(real, array[] real, array[] real, data array[] real, data array[] int):array[] real
```
@@ -140,7 +146,7 @@ derivatives with respect to time of the state,
For more fine-grained control of the ODE solvers, these parameters can
also be provided:
-
+
* `data` *`rel_tol`*: relative tolerance for the ODE solver, type
`real`, data only,
@@ -159,7 +165,7 @@ with values consisting of solutions at the specified times.
The sizes must match, and in particular, the following groups are of
the same size:
-
+
* state variables passed into the system function, derivatives
returned by the system function, initial state passed into the
solver, and rows of the return value of the solver,
@@ -169,3 +175,65 @@ solver,
* parameters, real data and integer data passed to the solver will
be passed to the system function
+
+
+## Exponentiated quadratic covariance functions {#cov_exp_quad}
+
+These covariance functions have been replaced by those described in:
+
+```{r results='asis', echo=FALSE}
+if (knitr::is_html_output()) {
+ cat(' * Gaussian Process Covariance Functions\n')
+}
+```
+
+With magnitude $\alpha$ and length scale $l$, the exponentiated quadratic kernel is:
+
+$$
+k(x_i, x_j) = \alpha^2 \exp \left(-\dfrac{1}{2\rho^2} \sum_{d=1}^D (x_{i,d} - x_{j,d})^2 \right)
+$$
+
+
+\index{{\tt \bfseries cov\_exp\_quad }!{\tt (row\_vectors x, real alpha, real rho): matrix}|hyperpage}
+
+`matrix` **`cov_exp_quad`**`(row_vectors x, real alpha, real rho)`
\newline
+The covariance matrix with an exponentiated quadratic kernel of x.
+`r since("2.16")`
+
+
+\index{{\tt \bfseries cov\_exp\_quad }!{\tt (vectors x, real alpha, real rho): matrix}|hyperpage}
+
+`matrix` **`cov_exp_quad`**`(vectors x, real alpha, real rho)`
\newline
+The covariance matrix with an exponentiated quadratic kernel of x.
+`r since("2.16")`
+
+
+\index{{\tt \bfseries cov\_exp\_quad }!{\tt (array[] real x, real alpha, real rho): matrix}|hyperpage}
+
+`matrix` **`cov_exp_quad`**`(array[] real x, real alpha, real rho)`
\newline
+The covariance matrix with an exponentiated quadratic kernel of x.
+`r since("2.16")`
+
+
+\index{{\tt \bfseries cov\_exp\_quad }!{\tt (row\_vectors x1, row\_vectors x2, real alpha, real rho): matrix}|hyperpage}
+
+`matrix` **`cov_exp_quad`**`(row_vectors x1, row_vectors x2, real alpha, real rho)`
\newline
+The covariance matrix with an exponentiated quadratic kernel of x1 and
+x2.
+`r since("2.18")`
+
+
+\index{{\tt \bfseries cov\_exp\_quad }!{\tt (vectors x1, vectors x2, real alpha, real rho): matrix}|hyperpage}
+
+`matrix` **`cov_exp_quad`**`(vectors x1, vectors x2, real alpha, real rho)`
\newline
+The covariance matrix with an exponentiated quadratic kernel of x1 and
+x2.
+`r since("2.18")`
+
+
+\index{{\tt \bfseries cov\_exp\_quad }!{\tt (array[] real x1, array[] real x2, real alpha, real rho): matrix}|hyperpage}
+
+`matrix` **`cov_exp_quad`**`(array[] real x1, array[] real x2, real alpha, real rho)`
\newline
+The covariance matrix with an exponentiated quadratic kernel of x1 and
+x2.
+`r since("2.18")`
diff --git a/src/functions-reference/matrix_operations.Rmd b/src/functions-reference/matrix_operations.Rmd
index 84710e342..d60545cb3 100644
--- a/src/functions-reference/matrix_operations.Rmd
+++ b/src/functions-reference/matrix_operations.Rmd
@@ -1604,72 +1604,314 @@ The cumulative sum of v
The cumulative sum of rv
`r since("2.0")`
-## Covariance functions {#covariance}
-
-### Exponentiated quadratic covariance function
-
-The exponentiated quadratic kernel defines the covariance between
-$f(x_i)$ and $f(x_j)$ where $f\colon \mathbb{R}^D \mapsto \mathbb{R}$
-as a function of the squared Euclidian distance between $x_i \in
-\mathbb{R}^D$ and $x_j \in \mathbb{R}^D$: \[ \text{cov}(f(x_i),
-f(x_j)) = k(x_i, x_j) = \alpha^2 \exp \left( -
-\dfrac{1}{2\rho^2} \sum_{d=1}^D (x_{i,d} - x_{j,d})^2 \right) \] with
-$\alpha$ and $\rho$ constrained to be positive.
-
-There are two variants of the exponentiated quadratic covariance
-function in Stan. One builds a covariance matrix, $K \in \mathbb{R}^{N
-\times N}$ for $x_1, \dots, x_N$, where $K_{i,j} = k(x_i, x_j)$, which
-is necessarily symmetric and positive semidefinite by construction.
-There is a second variant of the exponentiated quadratic covariance
-function that builds a $K \in \mathbb{R}^{N \times M}$ covariance
-matrix for $x_1, \dots, x_N$ and $x^\prime_1, \dots, x^\prime_M$,
-where $x_i \in \mathbb{R}^D$ and $x^\prime_i \in \mathbb{R}^D$ and
-$K_{i,j} = k(x_i, x^\prime_j)$.
-
-
-\index{{\tt \bfseries cov\_exp\_quad }!{\tt (row\_vectors x, real alpha, real rho): matrix}|hyperpage}
-
-`matrix` **`cov_exp_quad`**`(row_vectors x, real alpha, real rho)`
\newline
-The covariance matrix with an exponentiated quadratic kernel of x.
-`r since("2.16")`
-
-
-\index{{\tt \bfseries cov\_exp\_quad }!{\tt (vectors x, real alpha, real rho): matrix}|hyperpage}
-
-`matrix` **`cov_exp_quad`**`(vectors x, real alpha, real rho)`
\newline
-The covariance matrix with an exponentiated quadratic kernel of x.
-`r since("2.16")`
-
-
-\index{{\tt \bfseries cov\_exp\_quad }!{\tt (array[] real x, real alpha, real rho): matrix}|hyperpage}
-
-`matrix` **`cov_exp_quad`**`(array[] real x, real alpha, real rho)`
\newline
-The covariance matrix with an exponentiated quadratic kernel of x.
-`r since("2.16")`
-
-
-\index{{\tt \bfseries cov\_exp\_quad }!{\tt (row\_vectors x1, row\_vectors x2, real alpha, real rho): matrix}|hyperpage}
-
-`matrix` **`cov_exp_quad`**`(row_vectors x1, row_vectors x2, real alpha, real rho)`
\newline
-The covariance matrix with an exponentiated quadratic kernel of x1 and
-x2.
-`r since("2.18")`
+## Gaussian Process Covariance Functions {#gaussian-process-covariance-functions}
-
-\index{{\tt \bfseries cov\_exp\_quad }!{\tt (vectors x1, vectors x2, real alpha, real rho): matrix}|hyperpage}
+The Gaussian process covariance functions compute the covariance between
+observations in an input data set or the cross-covariance between two input
+data sets.
-`matrix` **`cov_exp_quad`**`(vectors x1, vectors x2, real alpha, real rho)`
\newline
-The covariance matrix with an exponentiated quadratic kernel of x1 and
-x2.
-`r since("2.18")`
+For one dimensional GPs, the input data sets are arrays of scalars. The
+covariance matrix is given by $K_{ij} = k(x_i, x_j)$ (where $x_i$ is the
+$i^{th}$ element of the array $x$) and the cross-covariance is given by
+$K_{ij} = k(x_i, y_j)$.
-
-\index{{\tt \bfseries cov\_exp\_quad }!{\tt (array[] real x1, array[] real x2, real alpha, real rho): matrix}|hyperpage}
+For multi-dimensional GPs, the input data sets are arrays of vectors. The
+covariance matrix is given by $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$ (where
+$\mathbf{x}_i$ is the $i^{th}$ vector in the array $x$) and the cross-covariance
+is given by $K_{ij} = k(\mathbf{x}_i, \mathbf{y}_j)$.
-`matrix` **`cov_exp_quad`**`(array[] real x1, array[] real x2, real alpha, real rho)`
\newline
-The covariance matrix with an exponentiated quadratic kernel of x1 and
-x2.
-`r since("2.18")`
+### Exponentiated quadratic kernel
+
+With magnitude $\sigma$ and length scale $l$, the exponentiated quadratic kernel is:
+
+$$
+k(\mathbf{x}_i, \mathbf{x}_j) = \sigma^2 \exp \left( -\frac{|\mathbf{x}_i - \mathbf{x}_j|^2}{2l^2} \right)
+$$
+
+\index{{\tt \bfseries gp\_exp\_quad\_cov }!{\tt (array[] real x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exp_quad_cov`**`(array[] real x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with exponentiated quadratic kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exp\_quad\_cov }!{\tt (array[] real x1, array[] real x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exp_quad_cov`**`(array[] real x1, array[] real x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with exponentiated quadratic
+kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exp\_quad\_cov }!{\tt (vectors x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exp_quad_cov`**`(vectors x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with exponentiated quadratic kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exp\_quad\_cov }!{\tt (vectors x, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exp_quad_cov`**`(vectors x, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process covariance with exponentiated quadratic kernel in multiple
+dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exp\_quad\_cov }!{\tt (vectors x1, vectors x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exp_quad_cov`**`(vectors x1, vectors x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with exponentiated quadratic
+kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exp\_quad\_cov }!{\tt (vectors x1, vectors x2, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exp_quad_cov`**`(vectors x1, vectors x2, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with exponentiated quadratic
+kernel in multiple dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+### Dot product kernel
+
+With bias $\sigma_0$ the dot product kernel is:
+
+$$
+k(\mathbf{x}_i, \mathbf{x}_j) = \sigma_0^2 + \mathbf{x}_i^T \mathbf{x}_j
+$$
+
+\index{{\tt \bfseries gp\_dot\_prod\_cov }!{\tt (array[] real x, real sigma): matrix}|hyperpage}
+
+`matrix` **`gp_dot_prod_cov`**`(array[] real x, real sigma)`
\newline
+
+Gaussian process covariance with dot product kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_dot\_prod\_cov }!{\tt (array[] real x1, array[] real x2, real sigma): matrix}|hyperpage}
+
+`matrix` **`gp_dot_prod_cov`**`(array[] real x1, array[] real x2, real sigma)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with dot product
+kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_dot\_prod\_cov }!{\tt (vectors x, real sigma): matrix}|hyperpage}
+
+`matrix` **`gp_dot_prod_cov`**`(vectors x, real sigma)`
\newline
+
+Gaussian process covariance with dot product kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_dot\_prod\_cov }!{\tt (vectors x1, vectors x2, real sigma): matrix}|hyperpage}
+
+`matrix` **`gp_dot_prod_cov`**`(vectors x1, vectors x2, real sigma)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with dot product
+kernel in multiple dimensions.
+`r since("2.20")`
+
+### Exponential kernel
+
+With magnitude $\sigma$ and length scale $l$, the exponential kernel is:
+
+$$
+k(\mathbf{x}_i, \mathbf{x}_j) = \sigma^2 \exp \left( -\frac{|\mathbf{x}_i - \mathbf{x}_j|}{l} \right)
+$$
+
+\index{{\tt \bfseries gp\_exponential\_cov }!{\tt (array[] real x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exponential_cov`**`(array[] real x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with exponential kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exponential\_cov }!{\tt (array[] real x1, array[] real x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exponential_cov`**`(array[] real x1, array[] real x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with exponential
+kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exponential\_cov }!{\tt (vectors x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exponential_cov`**`(vectors x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with exponential kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exponential\_cov }!{\tt (vectors x, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exponential_cov`**`(vectors x, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process covariance with exponential kernel in multiple
+dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exponential\_cov }!{\tt (vectors x1, vectors x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exponential_cov`**`(vectors x1, vectors x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with exponential
+kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_exponential\_cov }!{\tt (vectors x1, vectors x2, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_exponential_cov`**`(vectors x1, vectors x2, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with exponential
+kernel in multiple dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+### Matern 3/2 kernel
+
+With magnitude $\sigma$ and length scale $l$, the Matern 3/2 kernel is:
+
+$$
+k(\mathbf{x}_i, \mathbf{x}_j) = \sigma^2 \left( 1 + \frac{\sqrt{3}|\mathbf{x}_i - \mathbf{x}_j|}{l} \right) \exp \left( -\frac{\sqrt{3}|\mathbf{x}_i - \mathbf{x}_j|}{l} \right)
+$$
+
+\index{{\tt \bfseries gp\_matern32\_cov }!{\tt (array[] real x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern32_cov`**`(array[] real x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with Matern 3/2 kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern32\_cov }!{\tt (array[] real x1, array[] real x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern32_cov`**`(array[] real x1, array[] real x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with Matern 3/2
+kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern32\_cov }!{\tt (vectors x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern32_cov`**`(vectors x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with Matern 3/2 kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern32\_cov }!{\tt (vectors x, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern32_cov`**`(vectors x, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process covariance with Matern 3/2 kernel in multiple
+dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern32\_cov }!{\tt (vectors x1, vectors x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern32_cov`**`(vectors x1, vectors x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with Matern 3/2
+kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern32\_cov }!{\tt (vectors x1, vectors x2, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern32_cov`**`(vectors x1, vectors x2, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with Matern 3/2
+kernel in multiple dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+### Matern 5/2 kernel
+
+With magnitude $\sigma$ and length scale $l$, the Matern 5/2 kernel is:
+
+$$
+k(\mathbf{x}_i, \mathbf{x}_j) = \sigma^2 \left( 1 + \frac{\sqrt{5}|\mathbf{x}_i - \mathbf{x}_j|}{l} + \frac{5 |\mathbf{x}_i - \mathbf{x}_j|^2}{3l^2} \right)
+ \exp \left( -\frac{\sqrt{5} |\mathbf{x}_i - \mathbf{x}_j|}{l} \right)
+$$
+
+\index{{\tt \bfseries gp\_matern52\_cov }!{\tt (array[] real x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern52_cov`**`(array[] real x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with Matern 5/2 kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern52\_cov }!{\tt (array[] real x1, array[] real x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern52_cov`**`(array[] real x1, array[] real x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with Matern 5/2
+kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern52\_cov }!{\tt (vectors x, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern52_cov`**`(vectors x, real sigma, real length_scale)`
\newline
+
+Gaussian process covariance with Matern 5/2 kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern52\_cov }!{\tt (vectors x, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern52_cov`**`(vectors x, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process covariance with Matern 5/2 kernel in multiple
+dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern52\_cov }!{\tt (vectors x1, vectors x2, real sigma, real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern52_cov`**`(vectors x1, vectors x2, real sigma, real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with Matern 5/2
+kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_matern52\_cov }!{\tt (vectors x1, vectors x2, real sigma, array[] real length\_scale): matrix}|hyperpage}
+
+`matrix` **`gp_matern52_cov`**`(vectors x1, vectors x2, real sigma, array[] real length_scale)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with Matern 5/2
+kernel in multiple dimensions with a length scale for each dimension.
+`r since("2.20")`
+
+### Periodic kernel
+
+With magnitude $\sigma$, length scale $l$, and period $p$, the periodic kernel is:
+
+$$
+k(\mathbf{x}_i, \mathbf{x}_j) = \sigma^2 \exp \left(-\frac{2 \sin^2 \left( \pi \frac{|\mathbf{x}_i - \mathbf{x}_j|}{p} \right) }{l^2} \right)
+$$
+
+\index{{\tt \bfseries gp\_periodic\_cov }!{\tt (array[] real x, real sigma, real length\_scale, real period): matrix}|hyperpage}
+
+`matrix` **`gp_periodic_cov`**`(array[] real x, real sigma, real length_scale, real period)`
\newline
+
+Gaussian process covariance with periodic kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_periodic\_cov }!{\tt (array[] real x1, array[] real x2, real sigma, real length\_scale, real period): matrix}|hyperpage}
+
+`matrix` **`gp_periodic_cov`**`(array[] real x1, array[] real x2, real sigma, real length_scale, real period)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with periodic
+kernel in one dimension.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_periodic\_cov }!{\tt (vectors x, real sigma, real length\_scale, real period): matrix}|hyperpage}
+
+`matrix` **`gp_periodic_cov`**`(vectors x, real sigma, real length_scale, real period)`
\newline
+
+Gaussian process covariance with periodic kernel in multiple dimensions.
+`r since("2.20")`
+
+\index{{\tt \bfseries gp\_periodic\_cov }!{\tt (vectors x1, vectors x2, real sigma, real length\_scale, real period): matrix}|hyperpage}
+
+`matrix` **`gp_periodic_cov`**`(vectors x1, vectors x2, real sigma, real length_scale, real period)`
\newline
+
+Gaussian process cross-covariance of `x1` and `x2` with periodic
+kernel in multiple dimensions with a length scale for each dimension.
+`r since("2.20")`
## Linear algebra functions and solvers