---
title: "Lecture 14" 
subtitle: "Covariance Functions"
date: "10/24/2018"
fontsize: 11pt
output: 
  beamer_presentation:
    theme: metropolis
    highlight: pygments
    fig_caption: false
    keep_tex: true
    latex_engine: xelatex
    includes:
      in_header: ../settings.tex
---

```{r setup, include=FALSE}
library(dplyr)
library(ggplot2)

#devtools::install_github("thomasp85/patchwork")
library(patchwork)

set.seed(20180308)
```

```{r config, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  fig.width=7,
  fig.height=4.5,
  out.width="\\textwidth",
  fig.align="center",
  echo=TRUE,
  warning=FALSE
)

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

source("../util.R")
```


# PM2.5 Example - Revisited

## FRN Data

Measured PM2.5 data from an EPA monitoring station in Columbia, NJ.

```{r echo=FALSE}
load("../Lec01/data/frn_example.Rdata")

pm25 = pm25 %>%
  mutate(Date = lubridate::mdy(Date)) %>%
  mutate(day  = (Date - lubridate::mdy("1/1/2007") + 1) %>% as.integer()) %>% 
  select(-POC) %>%
  setNames(., tolower(names(.)))

qlm = lm(pm25~day+I(day^2), data=pm25)

ggplot(pm25, aes(x=date, y=pm25)) +
  geom_line() +
  geom_point()
```


## Previous JAGS Model

\scriptoutput
```{r}
gp_exp_model = "model{
  y ~ dmnorm(mu, inverse(Sigma))

  for (i in 1:N) {
    mu[i] <- beta[1]+ beta[2] * x[i] + beta[3] * x[i]^2
  }
  
  for (i in 1:(N-1)) {
    for (j in (i+1):N) {
      Sigma[i,j] <- sigma2 * exp(- pow(l*d[i,j],2))
      Sigma[j,i] <- Sigma[i,j]
    }
  }

  for (k in 1:N) {
    Sigma[k,k] <- sigma2 + sigma2_w
  }

  for (i in 1:3) {
    beta[i] ~ dt(coef[i], 2.5, 1)
  }
  sigma2_w ~ dnorm(5, 1/25) T(0,)
  sigma2   ~ dnorm(12.5, 1/25) T(0,)
  l        ~ dt(0, 2.5, 1) T(0,) 
}"
```

```{r echo=FALSE}
load(file="../Lec13/gp_jags.Rdata")
load(file="../Lec13/gp_pred.Rdata")
```


## Posterior Predictive Distribution

```{r echo=FALSE}
ggplot(pm25) +
  geom_line(aes(x=day, y=pm25)) +
  geom_point(aes(x=day, y=pm25)) +
  geom_ribbon(data=y_pred_df, aes(ymin=post_lower, ymax=post_upper, x=x_pred), fill="red", alpha=0.1)+
  geom_line(data=y_pred_df, aes(x=x_pred, y=post_med), color='red', size=1)
```



## Revised JAGS Model

\scriptoutput
```{r}
gp_exp_model = "model{
  y ~ dmnorm(mu, inverse(Sigma))

  for (i in 1:N) {
    mu[i] <- beta[1]+ beta[2] * x[i] + beta[3] * x[i]^2
  }
  
  for (i in 1:(N-1)) {
    for (j in (i+1):N) {
      Sigma[i,j] <- sigma2 * exp(- pow(l*d[i,j],2))
      Sigma[j,i] <- Sigma[i,j]
    }
  }

  for (k in 1:N) {
    Sigma[k,k] <- sigma2 + sigma2_w
  }

  for (i in 1:3) {
    beta[i] ~ dt(coef[i], 2.5, 1)
  }
  sigma2_w ~ dnorm(5, 1/25) T(0,)
  sigma2   ~ dnorm(12.5, 1/25) T(0,)
  
  eff_range ~ dunif(3,90)
  l = sqrt(3) / eff_range
}"
```

```{r echo=FALSE}
if (file.exists("gp_fit.rds")) {
  gp_fit = readRDS(file="gp_fit.rds")
} else {
  m = rjags::jags.model(
    textConnection(gp_exp_model), 
    data = list(
      y = pm25$pm25,
      x = pm25$day,
      d = dist(pm25$day) %>% as.matrix(),
      N = nrow(pm25),
      coef = coef(qlm)
    ),
    quiet = TRUE
  )
  
  update(m, n.iter=10000)
  
  gp_fit = rjags::coda.samples(
    m, variable.names=c("beta", "sigma2", "l", "sigma2_w", "eff_range"),
    n.iter=10000, thin=10
  )
  saveRDS(gp_fit, file="gp_fit.rds")
}
```

## Posterior - Betas

```{r echo=FALSE}
betas = tidybayes::gather_draws(gp_fit, beta[i]) %>%
  ungroup() %>%
  mutate(.variable = paste0(.variable, "[",i,"]")) %>%
  select(-i)

betas %>%
  group_by(.variable) %>%
  slice(seq(1,n(),length.out=500)) %>%
  ggplot(aes(x=.iteration, y=.value, color=.variable)) +
  geom_line() +
  facet_grid(.variable~., scales = "free_y")
```


## Posterior - Covariance Parameters

```{r echo=FALSE}
params = tidybayes::gather_draws(gp_fit, sigma2, sigma2_w, l, eff_range)

params %>%
  slice(seq(1,n(),length.out=500)) %>%
  ggplot(aes(x=.iteration, y=.value, color=.variable)) +
  geom_line() +
  facet_wrap(~.variable, scales="free_y")
```

## Posterior - Covariance Parameters

```{r echo=FALSE}
params %>%
  slice(seq(1,n(),length.out=500)) %>%
  ggplot(aes(x=.value, fill=.variable)) +
  geom_density() +
  facet_wrap(~.variable, scales="free") +
  guides(fill=FALSE)
```

## Posterior Predictive Distribution {.t}

```{r echo=FALSE}
if (file.exists("gp_pred.rds")) {
  gp_pred = readRDS("gp_pred.rds")
} else {
  x = pm25$day
  y = pm25$pm25
  
  n_post_samp = 1000
  
  x_pred = 1:365 + rnorm(365, 0, 0.01)
  y_pred = matrix(NA, nrow=n_post_samp, ncol=length(x_pred))
  colnames(y_pred) = paste0("Y_pred[", round(x_pred,0), "]")
  
  dist_o  = fields::rdist(x)
  dist_p  = fields::rdist(x_pred)
  dist_op = fields::rdist(x, x_pred)
  dist_po = t(dist_op)
  
  p = progress_estimated(n_post_samp)
  
  for(i in 1:n_post_samp) {
    l        = filter(params, .iteration == i, .variable == 'l')        %>% pull(.value)
    sigma2   = filter(params, .iteration == i, .variable == "sigma2")   %>% pull(.value)
    sigma2_w = filter(params, .iteration == i, .variable == "sigma2_w") %>% pull(.value)
    beta0    = filter(betas,  .iteration == i, .variable == "beta[1]")  %>% pull(.value)
    beta1    = filter(betas,  .iteration == i, .variable == "beta[2]")  %>% pull(.value)
    beta2    = filter(betas,  .iteration == i, .variable == "beta[3]")  %>% pull(.value)
    
    #cat(l,sigma2, sigma2_w, beta0, beta1, beta2,"\n", sep=", ")
    
    mu = beta0 + beta1*x + beta2*x^2
    mu_pred = beta0 + beta1*x_pred + beta2*x_pred^2
    
    cov_o  = sq_exp_cov(dist_o,  sigma2 = sigma2, l = l) + nugget_cov(dist_o,  sigma2 = sigma2_w)
    cov_p  = sq_exp_cov(dist_p,  sigma2 = sigma2, l = l) + nugget_cov(dist_p,  sigma2 = sigma2_w)
    cov_op = sq_exp_cov(dist_op, sigma2 = sigma2, l = l) + nugget_cov(dist_op, sigma2 = sigma2_w)
    cov_po = sq_exp_cov(dist_po, sigma2 = sigma2, l = l) + nugget_cov(dist_po, sigma2 = sigma2_w)
    
    inv = solve(cov_o, cov_op)
    cond_cov = cov_p - cov_po %*% inv
    cond_mu  = mu_pred + t(inv) %*% (y - mu)
    
    y_pred[i,] = cond_mu + t(chol(cond_cov)) %*% matrix(rnorm(length(x_pred)), ncol=1)
    
    p$tick()$print()
  }
  
  gp_pred = y_pred %>% post_summary() %>% mutate(x_pred = round(x_pred,0))
  saveRDS(gp_pred, file="gp_pred.rds")
}
```

```{r echo=FALSE}
ggplot(pm25) +
  geom_line(aes(x=day, y=pm25)) +
  geom_point(aes(x=day, y=pm25)) +
  geom_ribbon(data=gp_pred, aes(ymin=post_lower, ymax=post_upper, x=x_pred), fill="red", alpha=0.1)+
  geom_line(data=gp_pred, aes(x=x_pred, y=post_mean), color='red', size=1)
```





# More Covariance Functions

## Nugget Covariance {.t}
\vspace{-5mm}
$$ Cov(y_{t_i}, y_{t_j}) = \sigma^2 {1}_{\{h=0\}} \text{   where } h = |t_i - t_j|$$

\vspace{5mm}

```{r echo=FALSE, fig.height=3}
n_draw = 2
h = 0:20

cov = data_frame(h = h) %>%
  mutate(C = nugget_cov(d=h, sigma2 = 1))

cov_ex = dist(h) %>% 
  as.matrix() %>% 
  nugget_cov(sigma2 = 1) %>%
  rmvnorm(n_draw, sigma=.) %>% 
  as.data.frame() %>%
  setNames(paste0("Draw ", 1:n_draw)) %>%
  mutate(x = h) %>%
  tidyr::gather(draw, y, -x)
  


p1 = ggplot(cov, aes(x=h, y=C)) +
    geom_point()
p2 = ggplot(cov_ex, aes(x=x, y=y, col=draw)) +
      geom_line()

p1 + p2
```


## (- / Power / Square) Exponential Covariance {.t}
\vspace{-5mm}
$$ Cov(y_{t_i}, y_{t_j}) = \sigma^2\exp\left(-(h\,l)^p\right) \text{   where } h = |t_i - t_j|$$

```{r echo=FALSE, fig.height=4.5, fig.align="center"}
n_draw = 2
h = seq(0,1, len=100)

l = 12

cov = data_frame(h = h) %>%
  mutate(
    "Exp" = exp_cov(d=h, sigma2 = 1, l=l),
    "Pow Exp\n(p=1.5)" = pow_exp_cov(d=h, sigma2 = 1, l=l, p=1.5),
    "Sq Exp" = sq_exp_cov(d=h, sigma2 = 1, l=l)
  ) %>%
  tidyr::gather("Cov", "C", -h)


cov_exp = dist(h) %>% 
    as.matrix() %>% 
    exp_cov(sigma2 = 1, l=l) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x) %>% 
    mutate(model = "Exp Cov")
cov_pow_exp = dist(h) %>% 
    as.matrix() %>% 
    pow_exp_cov(sigma2 = 1, l=l, p=1.5) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x) %>% 
    mutate(model = "Pow Exp Cov")
cov_sq_exp = dist(h) %>% 
    as.matrix() %>% 
    sq_exp_cov(sigma2 = 1, l=l) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x) %>% 
    mutate(model = "Sq Exp Cov")

  


p1 = ggplot(cov, aes(x=h, y=C, color=Cov)) +
  geom_line() + 
  labs(title="Covariance - l=12, sigma2=1")
p2 = ggplot(cov_exp, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  guides(color=FALSE) +
  labs(title="Exponential")
p3 = ggplot(cov_pow_exp, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(color=FALSE) +
  labs(title="Powered Exponential (p=1.5)")
p4 = ggplot(cov_sq_exp, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(color=FALSE) +
  labs(title="Square Exponential")

p1 + p2 + p3 + p4 + plot_layout(ncol=2)
```

## Matern Covariance {.t}
\vspace{-10mm}
$$ Cov(y_{t_i}, y_{t_j}) = \sigma^2 ~ \frac{2^{1-\nu}}{\Gamma(\nu)} ~ \left(\sqrt{2\nu}\, h \cdot l\right)^\nu ~ K_\nu\left(\sqrt{2\nu} \, h \cdot l\right) \text{   where } h = |t_i - t_j|$$

```{r echo=FALSE, fig.height=4.5, fig.align="center"}
n_draw = 2
h = seq(0, 6, len=200)

l = 2

cov = data_frame(h = seq(0,6, len=1000)) %>%
  mutate(
    "v=1/2" = matern_cov(d=h, sigma2 = 1, l=l, nu=1),
    "v=3/2" = matern_cov(d=h, sigma2 = 1, l=l, nu=2),
    "v=5/2" = matern_cov(d=h, sigma2 = 1, l=l, nu=3)
  ) %>%
  tidyr::gather("Cov", "C", -h)


cov_mat1 = dist(h) %>% 
    as.matrix() %>% 
    matern_cov(sigma2 = 1, l=l, nu=1/2) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_mat2 = dist(h) %>% 
    as.matrix() %>% 
    matern_cov(sigma2 = 1, l=l, nu=3/2) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_mat3 = dist(h) %>% 
    as.matrix() %>% 
    matern_cov(sigma2 = 1, l=l, nu=3/2) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)

  

p1 = ggplot(cov, aes(x=h, y=C, color=forcats::as_factor(Cov))) +
  geom_line() + 
  labs(title="Covariance - l=2, sigma2=1",color="")
p2 = ggplot(cov_mat1, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  guides(colour=FALSE) +
  labs(title="Matern - v=1/2")
p3 = ggplot(cov_mat2, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Matern - v=3/2")
p4 = ggplot(cov_mat3, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Matern - v=5/2")

p1 + p2 + p3 + p4 + plot_layout(ncol=2)
```

## Matern Covariance {.t}

* $K_\nu$ is the modified Bessel function of the second kind.

\vspace{1mm}

* A Gaussian process with Matérn covariance has sample functions that are $\lceil \nu -1\rceil$ times differentiable.

\vspace{1mm}

* When $\nu = 1/2 + p$ for $p \in \mathbb{N}^+$ then the Matern has a simplified form (product of an exponential and a polynomial of order $p$).

\vspace{1mm}

* When $\nu = 1/2$ the Matern is equivalent to the exponential covariance.

\vspace{1mm}    

* As $\nu \to \infty$ the Matern converges to the square exponential covariance.



## Rational Quadratic Covariance {.t}
\vspace{-5mm}
$$ Cov(y_{t_i}, y_{t_j}) = \sigma^2 \left(1 + \frac{h^2 \, l^2}{\alpha}\right)^{-\alpha} \text{   where } h = |t_i - t_j|$$

```{r echo=FALSE, fig.height=4.5, fig.align="center"}
n_draw = 2
h = seq(0,1, len=100)

l = 12

cov = data_frame(h = h) %>%
  mutate(
    "alpha=1"  = rquad_cov(d=h, sigma2 = 1, l=l, a=1),
    "alpha=3"  = rquad_cov(d=h, sigma2 = 1, l=l, a=3),
    "alpha=10" = rquad_cov(d=h, sigma2 = 1, l=l, a=10)
  ) %>%
  tidyr::gather("Cov", "C", -h)


cov_rq1 = dist(h) %>% 
    as.matrix() %>% 
    rquad_cov(sigma2 = 1, l=l, a=1) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_rq2 = dist(h) %>% 
    as.matrix() %>% 
    rquad_cov(sigma2 = 1, l=l, a=3) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_rq3 = dist(h) %>% 
    as.matrix() %>% 
    rquad_cov(sigma2 = 1, l=l, a=10) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)

  


p1 = ggplot(cov, aes(x=h, y=C, color=forcats::as_factor(Cov))) +
  geom_line() + 
  labs(title="Covariance - l=12, sigma2=1", y="", color="")
p2 = ggplot(cov_rq1, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  guides(colour=FALSE) +
  labs(title="Rational Quadratic - alpha=1")
p3 = ggplot(cov_rq2, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Rational Quadratic - alpha=3")
p4 = ggplot(cov_rq3, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Rational Quadratic - alpha=10")

p1 + p2 + p3 + p4 + plot_layout(ncol=2)
```





## Rational Quadratic Covariance {.t}

* is a scaled mixture of squared exponential covariance functions with different characteristic length-scales ($l$).

\vspace{1mm}

* As $\alpha \to \infty$ the rational quadratic converges to the square exponential covariance.

\vspace{1mm}

* Has sample functions that are infinitely differentiable for any value of $\alpha$ 


## Spherical Covariance {.t}
\vspace{-10mm}
\footnotesize
$$ Cov(y_{t_i}, y_{t_j}) = \begin{cases}
\sigma^2\left(1 - \frac{3}{2} h \cdot l + \frac{1}{2} (h \cdot l)^3)\right) & \text{if   } 0 < h < 1/l \\
0 & \text{otherwise}
\end{cases} \text{   where } h = |t_i - t_j|$$

```{r echo=FALSE, fig.height=4.5, fig.align="center"}
n_draw = 2
h = seq(0, 1.1, len=100)

l = 1

cov = data_frame(h = h) %>%
  mutate(
    "l=1"   = sphere_cov(d=h, sigma2 = 1, l=1),
    "l=3" = sphere_cov(d=h, sigma2 = 1, l=3),
    "l=10" = sphere_cov(d=h, sigma2 = 1, l=10)
  ) %>%
  tidyr::gather("Cov", "C", -h)


cov_sph1 = dist(h) %>% 
    as.matrix() %>% 
    sphere_cov(sigma2 = 1, l=1) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_sph2 = dist(h) %>% 
    as.matrix() %>% 
    sphere_cov(sigma2 = 1, l=3) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_sph3 = dist(h) %>% 
    as.matrix() %>% 
    sphere_cov(sigma2 = 1, l=10) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)

  
p1 = ggplot(cov, aes(x=h, y=C, color=forcats::as_factor(Cov))) +
  geom_line() + 
  labs(title="Covariance - sigma2=1", y="", color="")
p2 = ggplot(cov_sph1, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  guides(colour=FALSE) +
  labs(title="Spherical - l=1")
p3 = ggplot(cov_sph2, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Spherical - l=3")
p4 = ggplot(cov_sph3, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Spherical - l=10")

p1 + p2 + p3 + p4 + plot_layout(ncol=2)
```

## Periodic Covariance {.t}
\vspace{-10mm}
$$ Cov(y_{t_i}, y_{t_j}) = \sigma^2 \exp\left(-2\, l^2 \sin^2\left(\pi\frac{h}{p}\right)\right)  \text{   where } h = |t_i - t_j| $$

```{r echo=FALSE, fig.height=4.5, fig.align="center"}
n_draw = 2
h = seq(0, 6, len=200)

l = 2

cov = data_frame(h = seq(0,4, len=1000)) %>%
  mutate(
    "p=1" = periodic_cov(d=h, sigma2 = 1, l=l, p=1),
    "p=2" = periodic_cov(d=h, sigma2 = 1, l=l, p=2),
    "p=3" = periodic_cov(d=h, sigma2 = 1, l=l, p=3)
  ) %>%
  tidyr::gather("Cov", "C", -h)


cov_per1 = dist(h) %>% 
    as.matrix() %>% 
    periodic_cov(sigma2 = 1, l=l, p=1) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_per2 = dist(h) %>% 
    as.matrix() %>% 
    periodic_cov(sigma2 = 1, l=l, p=2) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)
cov_per3 = dist(h) %>% 
    as.matrix() %>% 
    periodic_cov(sigma2 = 1, l=l, p=3) %>%
    rmvnorm(n_draw, sigma=.) %>% 
    as.data.frame() %>%
    setNames(paste0("Draw ", 1:n_draw)) %>%
    mutate(x = h) %>%
    tidyr::gather(draw, y, -x)

  


p1 = ggplot(cov, aes(x=h, y=C, color=forcats::as_factor(Cov))) +
  geom_line() + 
  labs(title="Covariance - l=2, sigma2=1", y="")
p2 = ggplot(cov_per1, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  guides(colour=FALSE) +
  labs(title="Periodic - p=1")
p3 = ggplot(cov_per2, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Periodic - p=2")
p4 = ggplot(cov_per3, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE) +
  labs(title="Periodic - p=3")

p1 + p2 + p3 + p4 + plot_layout(ncol=2)
```


## Linear Covariance {.t}
\vspace{-5mm}
$$ Cov(y_{t_i}, y_{t_j}) = \sigma^2_b + \sigma^2_v~(t_i-c)(t_j-c)$$

```{r echo=FALSE, fig.height=4, fig.align="center", out.width="0.7\\textwidth"}
x = seq(0,1, length=21)

n_draw = 3
cov_lin = linear_cov(x, c=0) %>% 
  rmvnorm(n_draw, sigma=.) %>% 
  as.data.frame() %>%
  setNames(paste0("Draw ", 1:n_draw)) %>%
  mutate(x = x) %>%
  tidyr::gather(draw, y, -x)

ggplot(cov_lin, aes(x=x, y=y, color=draw)) +
  geom_line() +
  guides(colour=FALSE)
```

## Combining Covariances {.t}

If we definite two valid covariance functions, $Cov_a(y_{t_i}, y_{t_j})$ and $Cov_b(y_{t_i}, y_{t_j})$ then the following are also valid covariance functions,

$$
\begin{aligned}
Cov_a(y_{t_i}, y_{t_j}) + Cov_b(y_{t_i}, y_{t_j}) \\
Cov_a(y_{t_i}, y_{t_j}) \times Cov_b(y_{t_i}, y_{t_j})
\end{aligned}
$$


## Linear $\times$ Linear $\to$ Quadratic {.t}
\vspace{-5mm}
$$ Cov_a(y_{t_i}, y_{t_j}) = 1 + 2~(t_i \times t_j) $$
$$ Cov_b(y_{t_i}, y_{t_j}) = 2 + 1~(t_i \times t_j) $$

```{r echo=FALSE, fig.height=4, fig.align="center"}
x = seq(-2,2, length=100)

n_draw = 3
cov_a = linear_cov(x, c=0, sigma2_b = 1, sigma2_v = 2)
cov_b = linear_cov(x, c=0, sigma2_b = 2, sigma2_v = 1)

cov_ll = (cov_a * cov_b) %>% 
  rmvnorm(n_draw, sigma=.) %>% 
  as.data.frame() %>%
  setNames(paste0("Draw ", 1:n_draw)) %>%
  mutate(x = x) %>%
  tidyr::gather(draw, y, -x)

ggplot(cov_ll, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  labs(title="Cov_a * Cov_b") +
  guides(colour=FALSE)
```

## Linear $\times$ Periodic {.t}

\vspace{-5mm}
$$ Cov_a(y_{t_i}, y_{t_j}) = 1 + 1~(t_i \times t_j) $$
$$ Cov_b(y_{t_i}, y_{t_j}) = \exp\left(-2\, \sin^2\left(2\pi\,h\right)\right) $$

```{r echo=FALSE, fig.height=4, fig.align="center"}
n_draw=2
x = seq(0,3, length=100)
d = dist(x) %>% as.matrix()

cov_a = linear_cov(x, c=0, sigma2_b = 1, sigma2_v = 2)
cov_b = periodic_cov(d, sigma2 = 1, l = 1, p = 0.5)

cov_lp = (cov_a * cov_b) %>% 
  rmvnorm(n_draw, sigma=.) %>% 
  as.data.frame() %>%
  setNames(paste0("Draw ", 1:n_draw)) %>%
  mutate(x = x) %>%
  tidyr::gather(draw, y, -x)

ggplot(cov_lp, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  labs(title="Cov_a * Cov_b") +
  guides(colour=FALSE)
```

## Linear + Periodic {.t}

\vspace{-5mm}
$$ Cov_a(y_{t_i}, y_{t_j}) = 1 + 1~(t_i \times t_j) $$
$$ Cov_b(h = |t_i - t_j|) = \exp\left(-2\, \sin^2\left(2\pi\,h\right)\right) $$

```{r echo=FALSE, fig.height=4, fig.align="center"}
n_draw=2
x = seq(0,3, length=100)
d = dist(x) %>% as.matrix()

cov_a = linear_cov(x, c=0, sigma2_b = 1, sigma2_v = 2)
cov_b = periodic_cov(d, sigma2 = 1, l = 1, p = 0.5)

cov_lp = (cov_a + cov_b) %>% 
  rmvnorm(n_draw, sigma=.) %>% 
  as.data.frame() %>%
  setNames(paste0("Draw ", 1:n_draw)) %>%
  mutate(x = x) %>%
  tidyr::gather(draw, y, -x)

ggplot(cov_lp, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  labs(title="Cov_a + Cov_b")
```

## Sq Exp $\times$ Periodic $\to$ Locally Periodic {.t}

\vspace{-5mm}
$$ Cov_a(h = |t_i - t_j|) =\exp(-(1/3)h^2) $$
$$ Cov_b(h = |t_i - t_j|) = \exp\left(-2\, \sin^2\left(\pi\,h\right)\right) $$

```{r echo=FALSE, fig.height=4, fig.align="center"}
n_draw=2
x = seq(0,6, length=100)
d = dist(x) %>% as.matrix()

cov_a = sq_exp_cov(d, sigma2 = 1, l = 1/3)
cov_b = periodic_cov(d, sigma2 = 1, l = 1, p = 1)

cov_sp = (cov_a * cov_b) %>% 
  rmvnorm(n_draw, sigma=.) %>% 
  as.data.frame() %>%
  setNames(paste0("Draw ", 1:n_draw)) %>%
  mutate(x = x) %>%
  tidyr::gather(draw, y, -x)

ggplot(cov_sp, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  guides(colour=FALSE) +
  labs(title="Cov_a * Cov_b")
```


## Sq Exp (short) + Sq Exp (long) {.t}
\vspace{-5mm}
$$ Cov_a(h = |t_i - t_j|) = (1/4) \exp(-4\sqrt{3}h^2) $$
$$ Cov_b(h = |t_i - t_j|) = \exp(-(\sqrt{3}/2)h^2) $$

```{r echo=FALSE, fig.height=4, fig.align="center"}
n_draw=2
x = seq(0,10, length=1000)
d = dist(x) %>% as.matrix()

cov_a = sq_exp_cov(d, sigma2 = 1/4, l = 4*sqrt(3))
cov_b = sq_exp_cov(d, sigma2 = 1, l = sqrt(3)/4)

cov_ss = (cov_a + cov_b) %>% 
  rmvnorm(n_draw, sigma=.) %>% 
  as.data.frame() %>%
  setNames(paste0("Draw ", 1:n_draw)) %>%
  mutate(x = x) %>%
  tidyr::gather(draw, y, -x)

ggplot(cov_ss, aes(x=x, y=y, color=draw)) +
  geom_line() + 
  guides(colour=FALSE) +
  labs(title="Cov_a + Cov_b")
```


## Sq Exp (short) + Sq Exp (long) (Seen another way) {.t}

```{r echo=FALSE, fig.height=5, fig.align="center"}
n_draw=2
x = seq(0,10, length=1000)
d = dist(x) %>% as.matrix()

cov_a = sq_exp_cov(d, sigma2 = 1/4, l = 4*sqrt(3))
cov_b = sq_exp_cov(d, sigma2 = 1, l = sqrt(3)/4)

y_a = rmvnorm(n_draw, sigma=cov_a) %>% as.data.frame() %>% setNames(paste0("Draw ", 1:n_draw))
y_b = rmvnorm(n_draw, sigma=cov_b) %>% as.data.frame() %>% setNames(paste0("Draw ", 1:n_draw))
y_c = y_a + y_b

rbind(
  y_a %>% mutate(name="Cov_A (short)", x=x),
  y_b %>% mutate(name="Cov_B (long)", x=x),
  y_c %>% mutate(name="Cov_A + Cov_B", x=x)
) %>%
  tidyr::gather(draw, y, starts_with("Draw")) %>%
  ggplot(aes(x=x, y=y, color=draw)) +
    geom_line() +
    guides(colour=FALSE) +
    facet_grid(draw~forcats::as_factor(name))
```

# BDA3 example

## BDA3

\begin{center}
\includegraphics[width=0.7\textwidth]{figs/bda_cover.png} 

$~$

\vspace{-3mm}
\url{http://research.cs.aalto.fi/pml/software/gpstuff/demo_births.shtml}
\end{center}



## Births (one year)

\begin{columns}
\begin{column}{0.5\textwidth}
\begin{center}
\includegraphics[width=\textwidth]{figs/births_pic1.png}
\end{center}
\end{column}
\begin{column}{0.5\textwidth}

1. Smooth long term trend \\ (\textit{sq exp cov})

\vspace{2mm}

2. Seven day periodic trend with decay (\textit{periodic $\times$ sq exp cov})

\vspace{2mm}

3. Constant mean

\end{column}
\end{columns}


## Component Contributions{.t}

We can view our GP in the following ways,

$$ \symbf{y} \sim \mathcal{N}(\symbf{\mu},~ \symbf{\Sigma}_1 + \symbf{\Sigma}_2 + \sigma^2 \symbf{I}\,) $$

but with appropriate conditioning we can also think of $\symbf{y}$ as being the sum of multipe independent GPs

$$ \symbf{y} = \symbf{\mu} + w_1(\symbf{t}) + w_2(\symbf{t}) + w_3(\symbf{t})$$
where 
$$
\begin{aligned}
w_1(\symbf{t}) &\sim \mathcal{N}(0, \symbf{\Sigma}_1) \\
w_2(\symbf{t}) &\sim \mathcal{N}(0, \symbf{\Sigma}_2) \\
w_3(\symbf{t}) &\sim \mathcal{N}(0, \sigma^2 \symbf{I}\,)
\end{aligned}
$$

## Decomposition of Covariance Components {.t}

\footnotesize
$$ \begin{bmatrix} 
y \\
w_1 \\
w_2
\end{bmatrix} 
\sim \mathcal{N} \left(
\begin{bmatrix}
\symbf{\mu} \\ 0 \\ 0
\end{bmatrix},~
\begin{bmatrix}
\Sigma_1 + \Sigma_2 + \sigma^2 \symbf{I} & \Sigma_1 & \Sigma_2 \\
\Sigma_1 & \Sigma_1 & 0\\
\Sigma_2 & 0  & \Sigma_2 \\
\end{bmatrix}
\right)
$$

\normalsize


therefore

$$ w_1 ~|~ \symbf{y},\symbf{\mu},\symbf{\theta} \sim \mathcal{N}(\symbf{\mu}_{cond},~ \symbf{\Sigma}_{cond}) $$

$$ \symbf{\mu}_{cond} = 0 + \Sigma_1 ~ (\Sigma_1 + \Sigma_2 + \sigma^2 I)^{-1}(\symbf{y}-\symbf{\mu}) $$
$$ \symbf{\Sigma}_{cond} = \Sigma_1 - \Sigma_1 (\Sigma_1 + \Sigma_2 + \sigma^2 \symbf{I})^{-1} {\Sigma_1}^t $$

## Births (multiple years)

\begin{adjustwidth}{-10mm}{-10mm}
\begin{center}
\includegraphics[width=\paperwidth]{figs/births_pic2.png}
\end{center}
\end{adjustwidth}

\vspace{-3mm}

\footnotesize

1. slowly changing trend (*sq exp cov*)

2. small time scale correlating noise (*sq exp cov*)

3. 7 day periodical component capturing day of week effect (*periodic $\times$ sq exp cov*)

4. 365.25 day periodical component capturing day of year effect (*periodic $\times$ sq exp cov*)

6. component to take into account the special days and interaction with weekends (*linear cov*)

7. independent Gaussian noise (*nugget cov*)

8. constant mean


# Mauna Loa Example

## Atmospheric CO$_2$

```{r echo=FALSE, message=FALSE}
loa = data_frame(
  y = co2 %>% strip_attrs(),
  x = time(co2) %>% strip_attrs()
)

noaa = readr::read_csv("../data/noaa_mauna_loa.csv") %>%
  transmute(x = year+(month-1)/12, y=average) %>%
  filter(x > max(loa$x))

rbind(
  loa %>% mutate(Source="Scripps (co2 in R)"),
  noaa %>% mutate(Source="NOAA")
) %>%
ggplot(aes(x=x,y=y, color=Source)) + 
  geom_line() 
```

## GP Model {.t}

Based on Rasmussen 5.4.3 (we are using slightly different data and parameterization)

$$ \symbf{y} \sim \mathcal{N}(\symbf{\mu},~ \symbf{\Sigma}_1 + \symbf{\Sigma}_2 + \symbf{\Sigma}_3 + \symbf{\Sigma}_4 + \sigma^2 \mathit{I}\,) $$

$$\{\symbf{\mu}\}_i = \bar{y}$$


$$
\begin{aligned}
\{\symbf{\Sigma}_1\}_{ij} &= \sigma^2_1 \exp\left(-(l_1 \cdot d_{ij})^2\right) \\
\{\symbf{\Sigma}_2\}_{ij} &= \sigma^2_2 \exp\left(-(l_2 \cdot d_{ij})^2\right)\exp\left(-2 \, (l_3)^2  \sin^2(\pi \, d_{ij} / p)\right) \\
\{\symbf{\Sigma}_3\}_{ij} &= \sigma^2_3 \left(1+\frac{(l_4 \cdot d_{ij})^2}{\alpha}\right)^{-\alpha} \\
\{\symbf{\Sigma}_4\}_{ij} &= \sigma^2_4 \exp\left(-(l_5 \cdot d_{ij})^2\right)
\end{aligned}
$$

## JAGS Model

\scriptoutput
```{r echo=TRUE}
ml_model = "model{
  y ~ dmnorm(mu, inverse(Sigma))

  for (i in 1:(length(y)-1)) {
    for (j in (i+1):length(y)) {
      k1[i,j] <- sigma2[1] * exp(- pow(l[1] * d[i,j],2))
      k2[i,j] <- sigma2[2] * exp(- pow(l[2] * d[i,j],2) - 2 * pow(l[3] * sin(pi*d[i,j] / per), 2))
      k3[i,j] <- sigma2[3] * pow(1+pow(l[4] * d[i,j],2)/alpha, -alpha)
      k4[i,j] <- sigma2[4] * exp(- pow(l[5] * d[i,j],2))
      
      Sigma[i,j] <- k1[i,j] + k2[i,j] + k3[i,j] + k4[i,j]
      Sigma[j,i] <- Sigma[i,j]
    }
  }

  for (i in 1:length(y)) {
    Sigma[i,i] <- sigma2[1] + sigma2[2] + sigma2[3] + sigma2[4] + sigma2[5]
  }  

  for(i in 1:5){
    sigma2[i] ~ dt(0, 2.5, 1) T(0,)
    l[i] ~ dt(0, 2.5, 1) T(0,)
  }
  alpha ~ dt(0, 2.5, 1) T(0,)
}"
```


```{r echo=FALSE}
if (!file.exists("results/init_model.Rdata")) {
  m = rjags::jags.model(
    textConnection(ml_model), 
    data = list(
      mu = mean(loa$y) * rep(1, nrow(loa)),
      y = loa$y,
      #x = loa$x,
      d = dist(loa$x) %>% as.matrix(),
      per = 1,
      pi = pi
    ),
    inits = list(
      sigma2 = c(66, 2.4, 0.66, 0.18, 0.19)^2,
      l = 1/c(67, 90, 1.3, 1.2, 1.6/12),
      alpha = 0.78
    ),
    n.adapt=10000
  )
  save(m, file="results/init_model.Rdata")
} else {
  load("results/init_model.Rdata")
}

if (!file.exists("results/fit_model.Rdata")) {
  gp_coda = coda.samples(
    m, variable.names=c("sigma2", "l", "alpha"),
    n.iter=10000,
    thin=10
  )
  save(gp_coda, file="results/fit_model.Rdata")
} else {
  load("results/fit_model.Rdata")
}


sigma2 = get_coda_parameter(gp_coda, "sigma")
l      = get_coda_parameter(gp_coda, "^l")
alpha  = get_coda_parameter(gp_coda, "alpha")
```

## Diagnostics

```{r echo=FALSE}
r = tidybayes::gather_samples(gp_coda, sigma2[i], l[i], alpha) %>%
  ungroup() %>%
  mutate(term = fixterm(term, i))
  
r %>%  
  mutate(term = forcats::as_factor(term)) %>%
  group_by(term) %>%
  slice(seq(1,n(),length.out = 300)) %>%
  ggplot(aes(x=.iteration, y=estimate)) +
    geom_line() +
    facet_wrap(~term, scale="free_y", ncol=5)
```


## Fit Components

```{r echo=FALSE}
if (!file.exists("results/ml_gp_fit_comp.Rdata"))
{
  x_pred = loa$x
  y_fit_comp = lapply(1:4, function(i) matrix(NA, nrow=n_post_samp, ncol=length(x_pred)))
    
  mu = rep(mean(y), length(x))
  mu_pred = rep(0, length(x_pred))
  
  dist_o = rdist(x)
  dist_p = rdist(x_pred)
  dist_op = rdist(x, x_pred)
  
  p = progress_estimated(n_post_samp)
  for(i in 1:nrow(sigma2))
  {
    cov_o  = calc_cov(dist_o, i)$S
    inv = solve(cov_o)
    
    for(j in 1:4)
    {
      cov_name = paste0("S",j)
      cov_p  = calc_cov(dist_p, i)[[cov_name]]
      cov_op = calc_cov(dist_op, i)[[cov_name]]
      
      cond_cov = cov_p - t(cov_op) %*% inv %*% cov_op
      diag(cond_cov) = diag(cond_cov) + 1e-4
      cond_mu  = mu_pred + t(cov_op) %*% inv %*% (y - mu)
    
      y_fit_comp[[j]][i,] = cond_mu + t(chol(cond_cov)) %*% matrix(rnorm(length(x_pred)), ncol=1)
    }
    p$tick()$print()
  }
  
  save(y_fit_comp, file="results/ml_gp_fit_comp.Rdata")
} else {
  load(file="results/ml_gp_fit_comp.Rdata")
}
```

```{r echo=FALSE}
lapply(
  1:4, 
  function(i) 
    post_summary(y_fit_comp[[i]]) %>% 
      mutate(cov=paste0("Sigma_",i), x=loa[["x"]])
) %>% 
  bind_rows() %>%
  ggplot(aes(x=x, y=post_mean, col=cov, fill=cov)) +
    geom_ribbon(aes(ymin=post_lower, ymax=post_upper, color=NULL), alpha=0.2) +
    geom_line() +
    facet_wrap(~cov, scale="free_y")
```




## Forecasting

```{r echo=FALSE}
y = loa$y
x = loa$x

n_post_samp = 1000

calc_cov = function(d,i)
{
  S1 = sigma2[i,1] * exp(- (l[i,1] * d)^2)
  S2 = sigma2[i,2] * exp(- (l[i,2] * d)^2 - 2 * (l[i,3] * sin(pi*d))^2)
  S3 = sigma2[i,3] * (1+(l[i,4] * d)^2/alpha[i,1])^-alpha[i,1]
  S4 = sigma2[i,4] * exp(- (l[i,5] * d)^2)
  S5 = ifelse(abs(d)<1e-4, sigma2[i,5] + 1e-6, 0)
  S = S1 + S2 + S3 + S4 + S5
  
  list(S1=S1, S2=S2, S3=S3, S4=S4, S5=S5, S=S)
}
```

```{r echo=FALSE}
if (!file.exists("results/ml_gp_fore.Rdata"))
{
  x_pred = noaa$x
  y_fore = matrix(NA, nrow=n_post_samp, ncol=length(x_pred))
  colnames(y_fore) = paste0("Y_pred[", 1:length(x_pred), "]")
  
  mu = rep(mean(y), length(x))
  mu_pred = rep(mean(y), length(x_pred))
  
  dist_o = fields::rdist(x)
  dist_p = fields::rdist(x_pred)
  dist_op = fields::rdist(x, x_pred)
  
  for(i in 1:nrow(sigma2))
  {
    cov_o  = calc_cov(dist_o, i)$S
    cov_p  = calc_cov(dist_p, i)$S
    cov_op = calc_cov(dist_op, i)$S
    
    inv = solve(cov_o, cov_op)
    cond_cov = cov_p - t(cov_op) %*% inv
    diag(cond_cov) = diag(cond_cov) + 1e-5
    cond_mu  = mu_pred + t(inv) %*% (y - mu)
    
    y_fore[i,] = cond_mu + t(chol(cond_cov)) %*% matrix(rnorm(length(x_pred)), ncol=1)
  }
  
  save(y_fore, file="results/ml_gp_fore.Rdata")
} else {
  load(file="results/ml_gp_fore.Rdata")
}
```

```{r echo=FALSE}
summ_80 = post_summary(y_fore, ci_width = 0.8) %>% mutate(x = noaa[["x"]])
summ_95 = post_summary(y_fore, ci_width = 0.95) %>% mutate(x = noaa[["x"]])

ggplot(summ_95, aes(x=x, y=post_mean)) +
  geom_ribbon(data=summ_95, aes(ymin=post_lower, ymax=post_upper), fill="#D3D2D6") +
  geom_ribbon(data=summ_80, aes(ymin=post_lower, ymax=post_upper), fill="#A3A5C3") +
  geom_line(color="#1007FC") +
  geom_line(data=loa, aes(x=x,y=y), alpha=0.5) +
  geom_line(data=noaa, aes(x=x,y=y), alpha=0.5, col='red')
```

## Forecasting (zoom)

```{r echo=FALSE}
ggplot(summ_95, aes(x=x, y=post_mean)) +
  geom_ribbon(data=summ_95, aes(ymin=post_lower, ymax=post_upper), fill="#D3D2D6") +
  geom_ribbon(data=summ_80, aes(ymin=post_lower, ymax=post_upper), fill="#A3A5C3") +
  geom_line(color="#1007FC", size=0.8) +
  geom_line(data=loa, aes(x=x,y=y), alpha=0.5, size=0.8) +
  geom_line(data=noaa, aes(x=x,y=y), alpha=0.5, col='red', size=0.8) +
  xlim(1998, 2017) + ylim(360,410)
```

## Forecasting ARIMA (auto)


```{r echo=FALSE}
m5 = forecast::auto.arima(co2)
m5_fore = forecast::forecast(m5,h = nrow(noaa))

autoplot(m5_fore) +
  geom_line(data=noaa, aes(x=x,y=y), alpha=0.5, col='red', size=0.8)
```



## Forecasting RMSE

```{r echo=FALSE}
noaa_pred_arima = data_frame(
  y_hat = m5_fore$mean %>% strip_attrs(),
  x = time(m5_fore$mean) %>% strip_attrs()
) %>% 
  mutate(y = noaa[["y"]])

noaa_pred_gp = summ_95 %>%
  select(y_hat = post_mean, x) %>%
  mutate(y = noaa[["y"]])

filters = c(2003, 2008, 2013, 2018)

noaa_rmse = function(year, d) 
{
  d %>% 
    filter(x < year) %>% 
    {(.$y - .$y_hat)^2} %>% 
    mean() %>% 
    sqrt() %>% 
    round(3)
}

data_frame(
  "dates"        = c("Jan 1998 - Jan 2003", "Jan 1998 - Jan 2008", "Jan 1998 - Jan 2013", "Jan 1998 - Mar 2017"),
  "RMSE (arima)" = purrr::map_dbl(filters, noaa_rmse, d=noaa_pred_arima),
  "RMSE (gp)"    = purrr::map_dbl(filters, noaa_rmse, d=noaa_pred_gp)
) %>%
  knitr::kable()
```

 


## Forecasting Components

```{r echo=FALSE}
if (!file.exists("results/ml_gp_fore_comp.Rdata"))
{
  x_pred = noaa$x
  y_comp = lapply(1:4, function(i) matrix(NA, nrow=n_post_samp, ncol=length(x_pred)))
    
  mu = rep(mean(y), length(x))
  mu_pred = rep(0, length(x_pred))
  
  dist_o = rdist(x)
  dist_p = rdist(x_pred)
  dist_op = rdist(x, x_pred)
  
  p = progress_estimated(n_post_samp)
  for(i in 1:nrow(sigma2))
  {
    cov_o  = calc_cov(dist_o, i)$S
    inv = solve(cov_o)
    
    for(j in 1:4)
    {
      cov_name = paste0("S",j)
      cov_p  = calc_cov(dist_p, i)[[cov_name]]
      cov_op = calc_cov(dist_op, i)[[cov_name]]
      
      cond_cov = cov_p - t(cov_op) %*% inv %*% cov_op
      diag(cond_cov) = diag(cond_cov) + 1e-4
      cond_mu  = mu_pred + t(cov_op) %*% inv %*% (y - mu)
    
      y_comp[[j]][i,] = cond_mu + t(chol(cond_cov)) %*% matrix(rnorm(length(x_pred)), ncol=1)
    }
    p$tick()$print()
  }
  
  save(y_comp, file="results/ml_gp_fore_comp.Rdata")
} else {
  load(file="results/ml_gp_fore_comp.Rdata")
}
```

```{r echo=FALSE}
lapply(
  1:4, 
  function(i) 
    post_summary(y_comp[[i]]) %>% 
      mutate(cov=paste0("Sigma_",i), x=noaa[["x"]])
) %>% 
  bind_rows() %>%
  ggplot(aes(x=x, y=post_mean, col=cov, fill=cov)) +
    geom_ribbon(aes(ymin=post_lower, ymax=post_upper, color=NULL), alpha=0.2) +
    geom_line() +
    facet_wrap(~cov, scale="free_y")
```

