---
title: "Gaussian Process Models"
subtitle: "Lecture 14"
author: "Dr. Colin Rundel"
footer: "Sta 344 - Fall 2022"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
execute:
  echo: true
  warning: false
  collapse: true
---

```{r setup, include=FALSE}
library(tidyverse)
library(tsibble)

knitr::opts_chunk$set(
  fig.align = "center"
)

options(width=50)

ggplot2::theme_set(ggplot2::theme_bw())

source("util-crps.R")
source("fix_pred_draws.R")

set.seed(20221019)
```

# Multivariate Normal

## Multivariate Normal Distribution

For an $n$-dimension multivate normal distribution with covariance $\boldsymbol{\Sigma}$ (positive semidefinite) can be written as

$$
\underset{n \times 1}{\boldsymbol{y}} \sim N(\underset{n \times 1}{\boldsymbol{\mu}}, \; \underset{n \times n}{\boldsymbol{\Sigma}})
$$



$$
\begin{pmatrix}
y_1\\ \vdots\\ y_n
\end{pmatrix}
\sim N\left(
\begin{pmatrix}
\mu_1\\ \vdots\\ \mu_n
\end{pmatrix}, \,
\begin{pmatrix}
\rho_{11}\sigma_1\sigma_1 & \cdots & \rho_{1n}\sigma_1\sigma_n \\
\vdots & \ddots & \vdots \\
\rho_{n1}\sigma_n\sigma_1 & \cdots & \rho_{nn}\sigma_n\sigma_n \\
\end{pmatrix}
\right)
$$


## Density

For the $n$ dimensional multivate normal given on the last slide, its density is given by

$$
(2\pi)^{-n/2} \; \det(\boldsymbol{\Sigma})^{-1/2} \; \exp{\left(-\frac{1}{2} \underset{1 \times n}{(\boldsymbol{y}-\boldsymbol{\mu})'} \underset{n \times n}{\boldsymbol{\Sigma}^{-1}} \underset{n \times 1}{(\boldsymbol{y}-\boldsymbol{\mu})}\right)} 
$$

and its log density is given by

$$
-\frac{n}{2} \log 2\pi - \frac{1}{2} \log \det(\boldsymbol{\Sigma}) - -\frac{1}{2} \underset{1 \times n}{(\boldsymbol{y}-\boldsymbol{\mu})'} \underset{n \times n}{\boldsymbol{\Sigma}^{-1}} \underset{n \times 1}{(\boldsymbol{y}-\boldsymbol{\mu})}
$$


## Sampling

To generate draws from an $n$-dimensional multivate normal with mean $\underset{n \times 1}{\boldsymbol{\mu}}$ and covariance matrix $\underset{n \times n}{\boldsymbol{\Sigma}}$, 

. . .

* Find a matrix $\underset{n \times n}{\boldsymbol{A}}$ such that $\boldsymbol{\Sigma} = \boldsymbol{A}\,\boldsymbol{A}^t$

    * most often we use $\boldsymbol{A} = \text{Chol}(\boldsymbol{\Sigma})$ where $\boldsymbol{A}$ is a lower triangular matrix.

. . .

* Draw $n$ iid unit normals, $N(0,1)$, as $\underset{n \times 1}{\boldsymbol{z}}$ 

. . .

* Obtain multivariate normal draws using
$$ \underset{n \times 1}{\boldsymbol{y}} = \underset{n \times 1}{\boldsymbol{\mu}} + \underset{n \times n}{\boldsymbol{A}} \, \underset{n \times 1}{\boldsymbol{z}} $$


## Bivariate Examples

$$ \boldsymbol{\mu} = \begin{pmatrix}0 \\ 0\end{pmatrix} \qquad \boldsymbol{\Sigma} = \begin{pmatrix}1 & \rho \\ \rho & 1 \end{pmatrix}$$

```{r echo=FALSE}
map_dfr(
  c(0.9,0.7,0.5,0.1,-0.1,-0.5,-0.7,-0.9),
  function(rho) {
    mvtnorm::rmvnorm(1000, sigma = matrix(c(1,rho,rho,1),2)) %>% 
      as_tibble(.name_repair = ~ c("x","y")) %>%
      mutate(param = paste0("rho=",rho))
  }
) %>%
  ggplot(aes(x=x,y=y)) +
    geom_density_2d() +
    geom_point(alpha=0.1, size=0.5) +
    facet_wrap(~param, ncol=4)
```

## Marginal / conditional distributions



*Proposition* - For an $n$-dimensional multivate normal with mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$, any marginal or conditional distribution of the $y$'s will also be (multivariate) normal.

. . .

*Univariate marginal distribution*:
$$ y_i = N(\boldsymbol{\mu}_i,\,\boldsymbol{\Sigma}_{ii}) $$

. . .

*Bivariate marginal distribution*:
$$ \boldsymbol{y}_{ij} = N\left( \begin{pmatrix}\boldsymbol{\mu}_i \\ \boldsymbol{\mu}_j \end{pmatrix},\; \begin{pmatrix} \boldsymbol{\Sigma}_{ii} & \boldsymbol{\Sigma}_{ij} \\ \boldsymbol{\Sigma}_{ji} & \boldsymbol{\Sigma}_{jj} \end{pmatrix} \right) $$

## 

*$k$-dimensional marginal distribution*:

$$ 
\boldsymbol{y}_{i,\ldots,k} = 
  N\left( 
    \begin{pmatrix}\boldsymbol{\mu}_{i} \\ \vdots \\ \boldsymbol{\mu}_{k} \end{pmatrix},\; 
    \begin{pmatrix} 
      \boldsymbol{\Sigma}_{ii}  & \cdots & \boldsymbol{\Sigma}_{i k} \\ 
      \vdots           & \ddots & \vdots \\
      \boldsymbol{\Sigma}_{k i} & \cdots & \boldsymbol{\Sigma}_{k k} \end{pmatrix} 
  \right) 
$$

::: {.aside}
$i,\ldots,k$ do not need to be ordered in any particular way
:::

## Conditional Distributions
 
If we partition the $n$-dimensions into two pieces such that $\boldsymbol{y} = (\boldsymbol{y}_1,\, \boldsymbol{y}_2)^t$ then

::: {.small}
$$
\underset{n \times 1}{\boldsymbol{y}} \sim N\left(
  \begin{pmatrix}
    \boldsymbol{\mu}_1 \\ 
    \boldsymbol{\mu}_2
  \end{pmatrix},\, 
  \begin{pmatrix} 
    \boldsymbol{\Sigma}_{11} & \boldsymbol{\Sigma}_{12} \\ 
    \boldsymbol{\Sigma}_{21} & \boldsymbol{\Sigma}_{22} 
  \end{pmatrix}
\right)
$$
$$ 
\begin{aligned}
\underset{k \times 1}{~~\boldsymbol{y}_1~~} &\sim N(\underset{k \times 1}{~~~\boldsymbol{\mu}_1~~~},\, \underset{k \times k}{~~~\boldsymbol{\Sigma}_{11}~~~}) \\ 
\underset{n-k \times 1}{\boldsymbol{y}_2} &\sim N(\underset{n-k \times 1}{\boldsymbol{\mu}_2},\, \underset{n-k \times n-k}{\boldsymbol{\Sigma}_{22}})
\end{aligned} 
$$
:::

. . .

then the conditional distributions are given by 

::: {.small}
$$
\begin{aligned}
\boldsymbol{y}_1 ~|~ \boldsymbol{y}_2 = \boldsymbol{a} ~&\sim N(\boldsymbol{\mu_1} + \boldsymbol{\Sigma_{12}} \, \boldsymbol{\Sigma_{22}}^{-1} \, (\boldsymbol{a} - \boldsymbol{\mu_2}),~ \boldsymbol{\Sigma_{11}}-\boldsymbol{\Sigma_{12}}\,\boldsymbol{\Sigma_{22}}^{-1} \, \boldsymbol{\Sigma_{21}}) \\
\\
\boldsymbol{y}_2 ~|~ \boldsymbol{y}_1 = \boldsymbol{b} ~&\sim N(\boldsymbol{\mu_2} + \boldsymbol{\Sigma_{21}} \, \boldsymbol{\Sigma_{11}}^{-1} \, (\boldsymbol{b} - \boldsymbol{\mu_1}),~ \boldsymbol{\Sigma_{22}}-\boldsymbol{\Sigma_{21}}\,\boldsymbol{\Sigma_{11}}^{-1} \, \boldsymbol{\Sigma_{12}})
\end{aligned}
$$
:::


## Gaussian Processes

From Shumway,

> A process, $\boldsymbol{y} = \{y(t) ~:~ t \in T\}$, is said to be a Gaussian process if all possible finite dimensional vectors $\boldsymbol{y} = (y_{t_1},y_{t_2},...,y_{t_n})^t$, for every collection of time points $t_1, t_2, \ldots , t_n$, and every positive integer $n$, have a multivariate normal distribution.

. . .

<br/>

So far we have only looked at examples of time series where $T$ is discete (and evenly spaces & contiguous), it turns out things get a lot more interesting when we explore the case where $T$ is defined on a *continuous* space  (e.g. $\mathbb{R}$ or some subset of $\mathbb{R}$).


# Gaussian Process Regression

## Parameterizing a Gaussian Process

Imagine we have a Gaussian process defined such that $\boldsymbol{y} = \{y(t) ~:~ t \in [0,1]\}$, 

. . .

* We now have an uncountably infinite set of possible $t$'s and $y(t)$s.

. . .

* We will only have a (small) finite number of observations $y(t_1), \ldots, y(t_n)$ with which to say something useful about this infinite dimensional process.

. . .

* The unconstrained covariance matrix for the observed data can have up to $n(n+1)/2$ unique values$^*$

. . .

* Necessary to make some simplifying assumptions:

    * Stationarity

    * Simple(r) parameterization of $\Sigma$


## Covariance Functions

More on these next week, but for now some common examples

. . .

*Exponential covariance*:
$$ \Sigma(y_{t},y_{t'}) = \sigma^2 \exp\big(-|t-t'| \; l\,\big) $$

. . .

*Squared exponential covariance (Gaussian)*:
$$ \Sigma(y_{t},y_{t'}) = \sigma^2 \exp\big(-(|t-t'| \; l\,)^2\big) $$

. . .

*Powered exponential covariance $\left(p \in (0,2]\right)$*:
$$ \Sigma(y_{t},y_{t'}) = \sigma^2 \exp\big(-(|t-t'| \; l\,)^p\big) $$

## Covariance Function - Correlation Decay

Letting $\sigma^2 = 1$ and trying different values of the length scale $l$,

```{r echo=FALSE, fig.height=4.5}
source("util-cov.R")
vals = expand_grid(
  d = seq(0, 1, length.out = 100),
  l = seq(1, 10, length.out = 10)
)

covs = rbind(
  mutate(vals, func="exp cov", corr = exp_cov(d,l=l)/exp_cov(0,l=l)),
  mutate(vals, func="sq exp cov", corr = sq_exp_cov(d,l=l)/sq_exp_cov(0,l=l))
)

ggplot(covs, aes(x=d, y=corr, color=as.factor(round(l,1)))) +
  geom_line() +
  facet_wrap(~func, ncol=2) +
  labs(color="l")
```

## Correlation Decay - AR(1)

Recall that for a stationary AR(1) process: $\gamma(h) = \sigma^2_w \phi^{|h|}$ and $\rho(h) = \phi^{|h|}$

. . .

we can draw a somewhat similar picture about the decay of $\rho$ as a function of distance.

```{r echo=FALSE}
expand_grid(
  phi = seq(0.1,0.9,0.2),
  h = 0:8
) %>%
  mutate(
    rho = phi^h,
    phi = as.factor(phi)
  ) %>%
  ggplot(aes(x=h, y=rho, color=phi)) +
    geom_point(size=2)

```


## Example

```{r echo=FALSE, fig.align="center"}
set.seed(1234)

y_f = function(x) log(x + 0.1) + sin(5*pi*x)

d_true = d = tibble(
  t = seq(0,1,length.out = 1000)
) %>%
  mutate(y = y_f(t) )

n = 20
d = tibble(
  t = seq(0,1,length.out = n) + rnorm(n, sd=0.01)
) %>%
  mutate(t = t - min(t)) %>%
  mutate(t = t / max(t)) %>%
  mutate(y = y_f(t) + rnorm(n,sd=0.2))

base = ggplot(d, aes(x=t, y=y)) +
  geom_point(size=3) +
  geom_line(data=d_true, color="blue", alpha=0.5, size=1)
base

save(d, d_true, y_f, base, file="lec14_ex.Rdata")
```

## Prediction

Our example has 20 observations which we would like to use as the basis for predicting $y(t)$ at other values of $t$ (say a regular sequence of values from 0 to 1).

. . .

<br/>

For now lets use a square exponential covariance with $\sigma^2 = 10$ and $l=15$ 

. . .

<br/>

We therefore want to sample from $\boldsymbol{y}_{pred} | \boldsymbol{y}_{obs}$.

$$\boldsymbol{y}_{pred} ~|~ \boldsymbol{y}_{obs} = \boldsymbol{y} ~\sim N(\boldsymbol{\Sigma}_{po} \, \boldsymbol{\Sigma}_{obs}^{-1} \, \boldsymbol{y},~ \boldsymbol{\Sigma}_{pred}-\boldsymbol{\Sigma}_{po}\,\boldsymbol{\Sigma}_{pred}^{-1} \, \boldsymbol{\Sigma}_{op})$$

```{r echo=FALSE}
d_pred = tibble(t = seq(0,1,length.out = 1000))
y_post = cond_pred(d_pred, d, cov=sq_exp_cov, sigma2 = 10, l = 15, reps=1000)

d_pred = cbind(
  d_pred,
  y_post[,1:5] %>% as.data.frame() %>% setNames(paste0("y",1:5)),
  post_summary(t(y_post))
)
```

## Draw 1

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), size=1)
```

## Draw 2

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=1.0, size=1)
```

## Draw 3

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=1.0, size=1)
```

## Draw 4

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y4), alpha=1.0, size=1)
```

## Draw 5

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y4), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y5), alpha=1.0, size=1)
```

## Many draws later

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y4), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y5), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=post_mean), size=1) +
  geom_ribbon(data=d_pred, fill='red', alpha=0.1, aes(ymin=post_lower, ymax=post_upper, y=post_mean))
```


## Exponential Covariance

Where $\sigma=10$, $l = \sqrt{15}$,

```{r echo=FALSE}
d_pred = tibble(t = seq(0,1,length.out = 1000))
y_post = cond_pred(d_pred, d, cov = exp_cov, sigma2 = 10, l = sqrt(15), reps=1000)

d_pred = cbind(
  d_pred,
  y_post[,1:5] %>% as.data.frame() %>% setNames(paste0("y",1:5)),
  post_summary(t(y_post))
)
```

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=1, size=0.5)
```

## Exponential Covariance - Draw 2

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=1, size=0.5)
```

## Exponential Covariance - Draw 3

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=1, size=0.5)
```

## Exponential Covariance - Draw 4

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y4), alpha=1, size=0.5)
```

## Exponential Covariance - Draw 5

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y4), alpha=0.2, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y5), alpha=1, size=0.5)
```


## Exponential Covariance - Variability

```{r echo=FALSE}
base + 
  geom_line(data=d_pred, color='red', aes(y=post_mean), size=1) +
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y4), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y5), alpha=0.4, size=0.5) +
  geom_ribbon(data=d_pred, fill='red', alpha=0.1, aes(ymin=post_lower, ymax=post_upper, y=post_mean))
```

## Powered Exponential Covariance ($p=1.5$)

```{r echo=FALSE}
d_pred = tibble(t = seq(0,1,length.out = 1000))

cov_15 = function(d, sigma2, l) pow_exp_cov(d,sigma2,l, p=1.5)

y_post = cond_pred(d_pred, d, cov = cov_15, sigma2 = 10, l = 15, reps=1000)

d_pred = cbind(
  d_pred,
  y_post[,1:5] %>% as.data.frame() %>% setNames(paste0("y",1:5)),
  post_summary(t(y_post))
)

base + 
  geom_line(data=d_pred, color='red', aes(y=post_mean), size=1) +
  geom_line(data=d_pred, color='red', aes(y=y1), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y2), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y3), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y4), alpha=0.4, size=0.5) +
  geom_line(data=d_pred, color='red', aes(y=y5), alpha=0.4, size=0.5) +
  geom_ribbon(data=d_pred, fill='red', alpha=0.1, aes(ymin=post_lower, ymax=post_upper, y=post_mean))
```

## Back to the square exponential

```{r echo=FALSE}
d_pred = tibble(t = seq(0,1,length.out = 100))

y_post = cond_pred(d_pred, d, cov = sq_exp_cov, sigma2 = 10, l = 15, reps=5000)

d_pred = cbind(
  d_pred,
  post_summary(t(y_post))
)

base + 
  geom_line(data=d_pred, color='red', aes(y=post_mean), size=1) +
  geom_ribbon(data=d_pred, fill='red', alpha=0.25, aes(ymin=post_lower, ymax=post_upper, y=post_mean))
```

## Changing the range ($l$)

```{r echo=FALSE}
d_pred = tibble(t = seq(0,1,length.out = 100))

d_pred_l = map_dfr(
  c(10,15,20,25),
  function(l) {
    cond_pred(d_pred, d, cov = sq_exp_cov, sigma2 = 10, l = l, reps=1000) %>% t() %>% 
      post_summary() %>%
      bind_cols(d_pred) %>%
      mutate(cov = paste0("Sq Exp Cov - sigma2=10, l=", l))
  }
) %>%
  mutate(cov = forcats::as_factor(cov))

base + 
  geom_line(data=d_pred_l, color='red', aes(y=post_mean), size=1) +
  geom_ribbon(data=d_pred_l, fill='red', alpha=0.25, aes(ymin=post_lower, ymax=post_upper, y=post_mean)) +
  facet_wrap(~cov)
```

## Effective Range

For the square exponential covariance
$$ \begin{aligned} 
Cov(d) &= \sigma^2 \exp\left(-(l \cdot d)^2\right) \\
Corr(d) &= \exp\left(-(l \cdot d)^2\right)
\end{aligned} $$

we would like to know, for a given value of $l$, beyond what distance apart must observations be to have a correlation less than $0.05$? 

. . .

$$
\begin{aligned}
\exp\left(-(l \cdot d)^2\right) &< 0.05 \\
-(l \cdot d)^2 &< \log 0.05 \\
l \cdot d &< \sqrt{3} \\
d &< \sqrt{3} / l
\end{aligned}
$$



## Changing the scale ($\sigma^2$)

```{r echo=FALSE}
d_pred = tibble(t = seq(0,1,length.out = 100))

d_pred_s = map2_dfr(
  c(5,5,5,15,15,15), c(10,15,20,10,15,20),
  function(s,l) {
    cond_pred(d_pred, d, cov = sq_exp_cov, sigma2 = s, l = l, reps=1000) %>% t() %>% 
      post_summary() %>%
      bind_cols(d_pred) %>%
      mutate(cov = paste0("Sq Exp Cov - sigma2=", s,", l=", l))
  }
) %>%
  mutate(cov = forcats::as_factor(cov))

base + 
  geom_line(data=d_pred_s, color='red', aes(y=post_mean), size=1) +
  geom_ribbon(data=d_pred_s, fill='red', alpha=0.25, aes(ymin=post_lower, ymax=post_upper, y=post_mean)) +
  facet_wrap(~cov)
```


## Fitting w/ BRMS

::: {.small}
```{r}
library(brms)
gp = brm(y ~ gp(t), data=d, cores=4, refresh=0)
```
:::

. . .

::: {.small}
```{r}
summary(gp)
```
:::

## Trace plots

```{r}
plot(gp)
```

## PP Checks

```{r}
pp_check(gp, ndraws = 100)
```

## Model predictions

```{r echo=FALSE}
ggplot(d, aes(x=t, y=y)) +
  tidybayes::stat_lineribbon(
    data = predicted_draws_fix(
      gp, newdata = tibble(t = seq(0,1,by=0.01))
    ),
    aes(y=.prediction),
    alpha=0.25
  ) +
  geom_line(
    data = d_true, color='blue'
  ) +
  geom_point()
```

## Forecasting

```{r echo=FALSE}
ggplot(d, aes(x=t, y=y)) +
  tidybayes::stat_lineribbon(
    data = predicted_draws_fix(
      gp, newdata = tibble(t = seq(0,1.5,by=0.01))
    ),
    aes(y=.prediction),
    alpha=0.25
  ) +
  geom_line(
    data = d_true, color='blue'
  ) +
  geom_point()
```

## Stan code

::: {.small}
```{r}
gp %>%
  brms::stancode()
```
:::