Consider Many Possible Functions

Rather than assuming some parametric model and then fitting parameters to it, in GPR we start by considering many possible functions that could occupy the x-y space. In our case, we can think of the functions as squiggly traces on a 2d plot. Before seeing any data, our possible functions could take just about any course or path in the space. Our first goal in GPR is to find a way to construct a bunch of random squiggly functions across a range of interest for x and over an assumed scale of y.

From Randomly Drawn Points to Traces

One way to construct a trace from left to right in 2d space would be to partition the x range up in to discrete intervals, draw a random y value at each x interval, and then connect the resulting points with lines. This would provide traces across the 2d space, but if the random draws were from a grid of y values or a uniform distribution in y, each draw’s position in y would be random and our trace would swing wildly up and down, displaying none of the cohesion or “smoothness” that we think of in functions. It would look like random noise.

n_uniforms <- 50
x_tbl <- tibble(x = seq(from = 0, to = 10, length = n_uniforms))

This function draws y values from a uniform [-3, 3] distribution at any specified x location. We then plot and take a look.

centered_unif_draws_fcn <- function(x_loc, n_draws){
  temp <- tibble(x = x_loc, y = runif(n_draws, -3, 3))
  return(temp)
}

unif_tbl <- x_tbl %>%
  mutate(y_draws = map(.x = x, .f = centered_unif_draws_fcn, n_draws = 600)) %>%
  select(-x) %>%
  unnest()

single_set <- x_tbl %>%
  mutate(y_draws = map(.x = x, .f = centered_unif_draws_fcn, n_draws = 1)) %>%
  select(-x) %>%
  unnest()

unif_tbl %>%
  ggplot(aes(x = x, y = y)) +
  geom_jitter(size = .1, alpha = .4, width = .02, color = "#2c3e50") +
  geom_point(data = single_set, aes(x = x, y = y), size = 2, color = "firebrick") +
  geom_line(data = single_set, aes(x = x, y = y), size = 1.5, color = "firebrick", alpha = .8) +
  labs(title = "Draws from Uniform Distributions",
       subtitle = "Positioned at x intervals",
       caption = "red line represents trace created by connecting a single draw from each distribtion")

Now here’s the amazing part - suppose that instead of drawing the y values from uniform distributions and creating noise, we use a normal (Gaussian) distribution instead. If we then asserted the Gaussians were actually part of a multi-variate normal distribution with a shared covariance matrix, we could control the nature of the draws such that points nearby each other (on the x axis) were highly correlated but points far from each other (on the x axis) were not highly correlated. It turns out that random points drawn and connected under these conditions yield a smooth trace.

A 2nd draw from the MVN produces another set of points that connect to form a 2nd random trace with a different path but similar smoothness and scale to the first. In this way, the stacking of normal curves in the x dimension allows random draws to produce random curves of specified scale and smoothness. We tailor the covariance matrix to give us the scale and smoothness to produce a range of various possible curves that could plausibly be the function that we want. These are our priors.

Uniform to Gaussian

This function draws y values from a standard normal distribution at any specified x location. We then plot and take a look.

centered_norm_draws_fcn <- function(x_loc){
  temp <- tibble(x = x_loc, y = rnorm(600, 0, 1))
  return(temp)
}

pretty_norm_tbl <- x_tbl %>%
  mutate(y_draws = map(.x = x, .f = centered_norm_draws_fcn)) %>%
  select(-x) %>%
  unnest() 

pp1 <- pretty_norm_tbl %>%
  ggplot(aes(x = x, y = y)) +
  geom_jitter(size = .1, alpha = .4, width = .02, color = "limegreen") 

pp2 <- pretty_norm_tbl %>%
  ggplot(aes(x = x, y = y)) +
  stat_halfeye(fill = "limegreen", alpha = .5)

pp1 / pp2 + plot_annotation(
  title = '(Multivariate) Normal Distributions Positioned at Intervals on the X Axis',
  subtitle = 'Each Normal Distribution has mean = 0 and sd = 1'
)

Calculate Distance Between any combination of points

The key idea to creating functions from random draws of a MVN is that points nearby each other are highly correlated while points farther away are less so. The distance() function (plgp library) calculates the distance of the partitioned x-range points from each other and outputs a matrix.

d_mat <- distance(x_tbl)

Interlude: A little nudge

We need to create a very small (negligible) numerical value to eventually add to the diagonal. This is a mathematical formality and things work fine without it in our simple case, but we add it here to be consistent with the reference materials. Don’t think too much about it.

diag_jitter <- sqrt(.Machine$double.eps) 
diag_jitter

## [1] 1.490116e-08

Get Covariance Matrix Sigma

We exponentiate the distances to give use the relationship between proximity of points to each other and their correlation. That’s all that is needed to make the covariance matrix!

sigma <- exp(-d_mat) + diag(diag_jitter, n_gaussians)

Draw Y-Values

Draw 1 value each from a bunch of normal distributions linked by the covariance matrix described above.

y_tbl <- tibble(y_draws = as.vector(rmvnorm(n = 1, mean = rep(0, length = nrow(x_tbl)), sigma = sigma))) 
y_tbl

## # A tibble: 50 x 1
##    y_draws
##      <dbl>
##  1   1.56 
##  2   1.49 
##  3   1.31 
##  4   1.08 
##  5   0.902
##  6   0.816
##  7   0.818
##  8   0.822
##  9   0.712
## 10   0.422
## # ... with 40 more rows

Combine the x’s and the y’s

Make a tibble to facilitate plotting.

gp_prior_tbl <- x_tbl %>%
  bind_cols(y_tbl)

gp_prior_tbl

## # A tibble: 50 x 2
##        x y_draws
##    <dbl>   <dbl>
##  1 0       1.56 
##  2 0.204   1.49 
##  3 0.408   1.31 
##  4 0.612   1.08 
##  5 0.816   0.902
##  6 1.02    0.816
##  7 1.22    0.818
##  8 1.43    0.822
##  9 1.63    0.712
## 10 1.84    0.422
## # ... with 40 more rows

Plot 1 Realization

By overlaying our randomly drawn points on the underlying Gaussians and connecting the dots, we see how the correlation matrix, based on distance in x, ensures that the points form a smooth but bumpy curve. The diagonal of the covariance matrix is 1, corresponding to ~95% of the green y values being between -2 and +2.

gp_prior_tbl %>%
  ggplot(aes(x = x, y = y_draws)) +
  geom_jitter(data = pretty_norm_tbl, aes(x = x, y = y), size = .1, alpha = .4, width = .02, color = "limegreen") +
  geom_point(size = 1, color = "black", alpha = .7) +
  geom_line(color = "black")+
  theme_classic() +
  labs(title = "Single Random Curve",
        subtitle = "Created Using Random Draw From MVN Distribution")

Define Function to Make a Random Trace

The above steps are converted to a function so that we can repeat our draws easily and create many traces.

generate_y_fcn <- function(sim_number){
  
x <- seq(from = 0, to = 10, length = n_gaussians)
d <- distance(x)
dj <- sqrt(.Machine$double.eps) 
sig <- exp(-d) + diag(dj, n_gaussians)
y <- tibble(x = x,
            y_draws = as.vector(rmvnorm(n = 1, mean = rep(0, length = length(x)), sigma = sig)),
            sim_number = sim_number)
return(y)
}

Map the Function Over a Setup Tbl to Create Multiple Traces

Generate and plot 12 instances of random curves over the background of Gaussians.

n_sim_traces <- 12

multi_trace_prior_tbl <- tibble(trace_num = seq(from = 1, to = n_sim_traces, by = 1) %>% as_factor()) %>%
  mutate(dat = map(.x = trace_num, .f = generate_y_fcn)) %>%
  unnest() 

multi_trace_prior_tbl %>% 
  ggplot(aes(x = x, y = y_draws)) +
  geom_jitter(data = pretty_norm_tbl, aes(x = x, y = y), size = .1, alpha = .9, width = .02, color = "grey") +
  geom_line(aes(color = trace_num), size = 1.3, alpha = .7) +
  geom_point(size = 1) +
  theme_classic() +
  theme(legend.position = "none") +
  scale_color_viridis_d(option = "Plasma") +
  labs(title = "Multiple Curves Representing Prior Possibilities",
       subtitle = "Created Using Random Draw From MVN Distribution")

Cool! We explored a clever way to create random functions from multi-normal draws and we specified a covariance matrix the controls the bumpiness and y-scale. Our priors are in place.

The Data Points

Generate a set of n=8 data point. In practice these would be supplied to us via an experiment or some other means.

n_dat <- 8
x_dat <- matrix(seq(0, 2*pi, length=n_dat), ncol=1)
y_dat <- sin(x_dat)

data_tbl <- tibble(x = x_dat, 
                   y = y_dat)

data_tbl %>%
  gt()

Calculating the Necessary Distances and Covariance Matrices

d_dat <- distance(x_dat) 
sigma_dat <- exp(-d_dat) + diag(diag_jitter, ncol(d_dat))
inv_sigma_dat <- solve(sigma_dat) #solve() finds inverse

x_pred_range <- matrix(seq(-0.5, 2*pi + 0.5, length=100), ncol=1)

d_x_pred_range <- distance(x_pred_range)
sigma_x_pred_range <- exp(-d_x_pred_range) + diag(diag_jitter, ncol(d_x_pred_range))

d_x_dat <- distance(x_pred_range, x_dat)
sigma_x_dat <- exp(-d_x_dat) 

Combine Per the Conditioning Equations Above

Bring it together by combining the above inputs into the predictive equations. The output will be trained means and (reduced) variances from which to draw new points to create new traces.

mu_pred <- sigma_x_dat %*% inv_sigma_dat %*% y_dat

sigma_pred <- sigma_x_pred_range - sigma_x_dat %*% inv_sigma_dat %*% t(sigma_x_dat)

Draw from the Updated Gaussians

y_draws <- rmvnorm(1, mu_pred, sigma_pred)

y_draws %>% as.vector()

##   [1] -0.11135178 -0.13414543 -0.14125646 -0.13384335 -0.11310097 -0.08032221
##   [7] -0.03869301  0.01008520  0.06426747  0.12167776  0.18125818  0.24235780
##  [13]  0.30473905  0.36878655  0.43483034  0.50286496  0.57276290  0.64341048
##  [19]  0.71394154  0.78186652  0.84577472  0.90297749  0.95168825  0.99017006
##  [25]  1.01800993  1.03518813  1.04214564  1.04042284  1.03171212  1.01776905
##  [31]  0.99957608  0.97924481  0.95714576  0.93327439  0.90739967  0.87885348
##  [37]  0.84624283  0.80916023  0.76606525  0.71690098  0.66145810  0.59996973
##  [43]  0.53349461  0.46304360  0.38937083  0.31415990  0.23807037  0.16243279
##  [49]  0.08706189  0.01247238 -0.06049590 -0.13256227 -0.20299983 -0.27120336
##  [55] -0.33690887 -0.39869342 -0.45653362 -0.50974602 -0.55842531 -0.60249705
##  [61] -0.64369848 -0.68291472 -0.72150414 -0.76086850 -0.80247216 -0.84637868
##  [67] -0.89223759 -0.93883876 -0.98369467 -1.02455925 -1.05798305 -1.08237703
##  [73] -1.09435116 -1.09289825 -1.07770612 -1.04977391 -1.00965956 -0.96046396
##  [79] -0.90481696 -0.84427854 -0.78205657 -0.71951229 -0.65718008 -0.59550566
##  [85] -0.53379511 -0.47180457 -0.40834782 -0.34347059 -0.27688415 -0.20952401
##  [91] -0.14200053 -0.07558342 -0.01236511  0.04618189  0.09857431  0.14360497
##  [97]  0.18013374  0.20735434  0.22425309  0.23046425

Plot the Draws

Plot the single trace drawn from posterior.

tibble(x = x_pred_range,
       y = as.vector(y_draws)) %>%
  ggplot(aes(x = x, y = y)) +
  geom_point(color = "purple") +
  geom_line(color = "purple") +
  labs(title = "Single Trace Generated from Draws of Posterior",
       subtitle = "Gaussians conditioned on n=8 data points")

Function for Many Posterior Draws

Now that we’ve shown how to create a set of draws to construct a posterior trace, let’s make it a function so that we can plot many traces. The function below is just all the same steps of the last section but combined into a tibble as the last step. This let’s us map the function easily as we expand the simulation.

generate_ypred_fcn <- function(sim_number){

#distance and sigma
d_dat <- distance(x_dat) 
sigma_dat <- exp(-d_dat) + diag(diag_jitter, ncol(d_dat))

#x range
x_pred_range <- matrix(seq(-0.5, 2*pi + 0.5, length=50), ncol=1)

#x point distance and sigma
d_x_pred_range <- distance(x_pred_range)
sigma_x_pred_range <- exp(-d_x_pred_range) + diag(diag_jitter, ncol(d_x_pred_range))

#distance and sigma between x points and x data points
d_x_dat <- distance(x_pred_range, x_dat)
sigma_x_dat <- exp(-d_x_dat) 

#inverse sigma for distance between x data points
inv_sigma_dat <- solve(sigma_dat)

#mu and sigma pred
mu_pred <- sigma_x_dat %*% inv_sigma_dat %*% y_dat
sigma_pred <- sigma_x_pred_range - sigma_x_dat %*% inv_sigma_dat %*% t(sigma_x_dat)

#y preds at x pred range points
y_draws <- rmvnorm(1, mu_pred, sigma_pred)

y <- tibble(x = x_pred_range,
            mu_y_pred = as.vector(mu_pred),
            y_draws = as.vector(rmvnorm(n = 1, mean = mu_pred, sigma = sigma_pred)),
            sim_number = sim_number)
return(y)
}

Generate Many Draws for Posterior Traces

Simply specify the number of desired traces and use the newly created function to start drawing that many from the MVN.

n_sim_traces <- 500

multi_trace_post_tbl <- tibble(trace_num = seq(from = 1, to = n_sim_traces, by = 1) %>% as_factor()) %>%
  mutate(dat = map(.x = trace_num, .f = generate_ypred_fcn)) %>%
  unnest() %>%
  select(-sim_number)

multi_trace_post_tbl %>%
  gt_preview()

Visualize

Structuring the data this way allows for plotting many traces of the posterior draws. In the figure below, the black line represents predictions, the red points are the training data, and the “bubbles” represent uncertainty. Notice how the uncertainty collapses near the observed data - just like we would want for FEA!

multi_trace_post_tbl %>% 
  ggplot(aes(x = x, y = y_draws)) +
  geom_line(aes(color = trace_num)) +
  geom_line(aes(x = x, y = mu_y_pred), size = 2) +
  geom_point(size = .5) +
  geom_point(dat = data_tbl, aes(x = x, y = y), size = 4, color = "firebrick") +
  theme(legend.position = "none") +
  scale_color_viridis_d() +
  labs(title = "Credible Posterior Traces",
       subtitle = "Gaussians trained on n=8 data points",
       caption = "Red points are observed data\nBlack line is predictions")

More Visualize

Just for fun, let’s just observe how the Gaussians collapse from prior to posterior, based on the data:

x_pred_range_new <- matrix(seq(-0.5, 2*pi + 0.5, length=50), ncol=1)

pretty_norm_tbl_2 <- tibble(x = x_pred_range_new) %>%
  mutate(y_draws = map(.x = x, .f = centered_norm_draws_fcn)) %>%
  select(-x) %>%
  unnest() 

prior_gaussians <- pretty_norm_tbl_2 %>% 
  ggplot(aes(x = x, y = y)) +
  geom_jitter(size = .1, alpha = .4, width = .01, color = "limegreen") 


post_gaussians <- multi_trace_post_tbl %>% 
  ggplot(aes(x = x, y = y_draws)) +
  geom_jitter(size = .1, alpha = .4, width = .01, color = "firebrick") +
  geom_point(dat = data_tbl, aes(x = x, y = y), size = 2, color = "black") 

prior_gaussians / post_gaussians +
  plot_annotation(title = "Gaussian Process Regression",
                  subtitle = "Image represents stacking of multivariate normal dimensions across an X-range of interest",
                  caption = "Upper Panel: Prior\nLower Panel: Posterior")