---
title: "Week 8: Activity"
format: pdf
editor: source
editor_options: 
  chunk_output_type: console
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.height = 6, fig.width = 6)
library(tidyverse) 
library(knitr)
library(mnormt)
library(plgp)
library(reshape2)

```


### Last Week's Recap

- Exam 1

- Linear Algebra 
- Linear Model Overview
- Simulating Data in R
- Fitting Linear Models in R
- Multivariate Normal distribution
- Partitioned Matrices and Conditional Multivariate normal distribution

### Video Lectures

- `rstan` in R for Bayesian inference

### This week

- Gaussian Process Intro
- Bayesian inference with `stan`
- Correlation functions

---


#### Visual Example

Let $n_1 = 1$ and $n_2 = 1$, then

$$\begin{bmatrix}y_1 \\ y_2 \end{bmatrix} \sim N \left( \begin{bmatrix}\mu_1 \\ \mu_2 \end{bmatrix}   ,
\begin{bmatrix} \sigma^2_{1} & \sigma_{12} \\ \sigma_{12} & \sigma_{2}^2 \end{bmatrix}\right)$$

and

$$y_1|y_2 \sim N \left( \mu_1 + \sigma_{12} (\sigma_{2}^{2})^{-1}\left(y_2 - \mu_2 \right), \sigma_{1}^2 - \sigma_{12} (\sigma_{2}^{2})^{-1} \sigma_{21} \right)$$
\vfill

```{r, echo = F}
#| fig-height: 4

mu1 <- 0
mu2 <- 0
sigmasq1 <- sigmasq2 <- 1

dat_seq <- seq(-4,4, by = .01)
n_seq <- length(dat_seq)
tibble(group = rep(c('y2 = -2; sigma12 = 0',
                     'y2 = -2; sigma12 = .2',
                     'y2 = -2; sigma12 = .8'), each = n_seq), 
dens = c(dnorm(dat_seq, mu1 + 0*(1/sigmasq2)*(-2 - mu2),
               sqrt(sigmasq1 - 0 * (1/sigmasq2)*0 )), 
         dnorm(dat_seq, mu1 + .2*(1/sigmasq2)*(-2 - mu2),
               sqrt(sigmasq1 - .2 * (1/sigmasq2)*.2 )), 
         dnorm(dat_seq, mu1 + .8*(1/sigmasq2)*(-2 - mu2),
               sqrt(sigmasq1 - .8 * (1/sigmasq2)*.8 ))), 
y = rep(dat_seq, 3)) %>% 
  ggplot(aes(x=y, y = dens, group = group, color = group)) +
  geom_line() + theme_bw()
```

\vfill

__Q:__ Calculate and write out the actual distributions for $y_1$ in these three settings.



\vfill

One last note, the marginal distributions for any partition $\underline{y_1}$ are quite simple.

$$\underline{y_1} \sim N \left( X_1\beta, \Sigma_{11} \right)$$
or just

$$y_1 \sim N \left( X_1\beta, \sigma^2_{1} \right)$$
if $y_1$ is scalar.


#### GP Overview

Now let's extend this idea to a Gaussian Process (GP). There are two fundamental ideas to a GP.

\vfill

1. Any finite set of realizations (say $\underline{y_2}$) has a multivariate normal distribution.

\vfill

2. Conditional on a set of realizations, all other locations (say $\underline{y_1}$) have a conditional normal distribution characterized by the mean, and most importantly the covariance function. Note the dimension of $\underline{y_1}$ can actually be infinite, such as defined on the real line.

\vfill

The big question is how to we estimate $\Sigma_{12}$? How many parameters are necessary for this distribution?

\vfill




#### Correlation function

Initially, let's consider correlation as a function of distance, in one dimension or on a line.


As a starting point, consider a variant of what is known as the exponential covariance function - we used this earlier. First define $d$ as the Euclidean distance between $x_1$ and $x_2$, such that $d = \sqrt{(x_i - x_j)^2}$ 

$$\rho_{i,j} = \exp \left(- d \right)$$

\vfill

Create a figure that shows the exponential correlation as a function of distance between the two points.

\vfill

Using a correlation function can reduce the number of unknown parameters in a covariance matrix. In an unrestricted case, $\Sigma$ has $n \choose 2$ + $n$ unknown parameters. However, using a correlation function can reduce the number of unknown parameters substantially, generally less than 4.


#### Realizations of a Gaussian Process

Recall that a process implies an infinite dimensional object. So we can generate a line rather than a discrete set of points. (While in practice the line will in fact be generated with a discrete set of points and then connected.)

\vfill

For this scenario we will assume a zero-mean GP, with covariance equal to the correlation function using $\rho_{i,j} = \exp \left(- d \right)$


\vfill

```{r}
#| echo: false
set.seed(02252025)
x <- seq(0, 10, by = .1)

dist_mat <- as.matrix(dist(x, diag = T, upper = T))
Sigma <- exp(-dist_mat)

y <- rmnorm(n =1, mean = 0 , varcov = Sigma)

tibble(y = y, x = x) |>
  ggplot(aes(y=y, x=x)) + 
  theme_bw() +
  geom_line() + 
  geom_hline(yintercept = 0, color = 'grey40') +
  ggtitle('Random realization of a GP') 
```

Overlay a few realizations of a Gaussian process on the same curve.




#### Connecting a GP to conditional normal

Now consider a discrete set of points, say $\underline{y_2}$, how can we estimate the response for the remainder of the values in the interval [0,10].


We can connect the dots (with uncertainty) using:

$$\underline{y_1}|\underline{y_2} \sim N \left( X_1\beta + \Sigma_{12} \Sigma_{22}^{-1}\left(\underline{y_2} - X_2\beta \right), \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21} \right)$$
\vfill

__Create a figure that shows the data points, conditional mean and uncertainty__


### GP Regression

Now rather than specifying a zero-mean GP, let the mean be $X \underline{\beta}.$ Create this figure.


---
#### Correlation function: more details

Recall the variant of the exponential covariance function that we have previously seen. Where $d$ as the Euclidean distance between $x_1$ and $x_2$, such that $d = \sqrt{(x_i - x_j)^2}$

$$\rho_{i,j} = \exp \left(- d \right)$$

\vfill

Recall that we can view the exponential correlation as a function of distance between the two points.


Now let's consider a more general framework for covariance where 

$$\sigma_{i,j} = \sigma^2 \exp \left(- d_{ij} /\phi \right)$$

Now we have introduced two new parameters into this function. What do you suppose that they do?

- $\sigma^2$: 

\vfill

- $\phi$: 

\vfill

Modify your previous code do adjust $\phi$ and $\sigma^2$ and explore how they differ.\vfill




We will soon talk about a more broad set of correlation functions and another parameter that provides flexibility so that predictions do not have to directly through observed points.

\vfill

### Geostatistical Data Exercise

At last, we will look at simulated 2-d "spatial" data.


##### 1. Create Sampling Locations


##### 2. Calculate Distances


##### 3. Define Covariance Function and Set Parameters

##### 4. Sample a realization of the process

- This requires a distributional assumption, we will use the Gaussian distribution

- Start with a coarse grid and them move to a finer grid

How does the spatial process change with:

- another draw with same parameters?
- a different value of $\phi$
- a different value of $\sigma^2$

