16 Day 16 (March 23)

16.1 Announcements

  • Assignment 7 is posted and due Sunday (March 26)
    • Derived quantities
  • Read Ch. 17 on spatial models

16.2 Gaussian process

  • A Gaussian process is a probability distribution over functions
    • If the function is observed at a finite number of points or “locations,” then the vector of values follows a multivariate normal distribution.

16.2.1 Multivariate normal distribution

  • Multivariate normal distribution

    • \(\boldsymbol{\eta}\sim\text{N}(\mathbf{0},\sigma^{2}\mathbf{R})\)
    • Definitions
      • Correlation matrix – A positive semi-definite matrix whose elements are the correlation between observations
      • Correlation function – A function that describes the correlation between observations
    • Example correlation matrices
      • Compound symmetry \[\mathbf{R}=\left[\begin{array}{cccccc} 1 & 1 & 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 1 & 0 & 0\\ 0 & 0 & 1 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 0 & 1 & 1 \end{array}\right]\]
      • AR(1) \[\mathbf{R(\phi})=\left[\begin{array}{ccccc} 1 & \phi^{1} & \phi^{2} & \cdots & \phi^{n-1}\\ \phi^{1} & 1 & \phi^{1} & \cdots & \phi^{n-2}\\ \phi^{2} & \phi^{1} & 1 & \vdots & \phi^{n-3}\\ \vdots & \vdots & \vdots & \ddots & \ddots\\ \phi^{n-1} & \phi^{n-2} & \phi^{n-3} & \ddots & 1 \end{array}\right]\:\]

  • Example simulating from \(\boldsymbol{\eta}\sim\text{N}(\mathbf{0},\sigma^{2}\mathbf{R})\) in R

    n <- 200
x <- 1:n
I <- diag(1, n)
sigma2 <- 1

library(MASS)
set.seed(2034)
eta <- mvrnorm(1, rep(0, n), sigma2 * I)
cor(eta[1:(n - 1)], eta[2:n])
    ## [1] -0.06408623
    acf(eta)

    par(mfrow = c(2, 1))
par(mar = c(0, 2, 0, 0), oma = c(5.1, 3.1, 2.1, 2.1))
plot(x, eta, typ = "l", xlab = "", xaxt = "n")
abline(a = 0, b = 0)

n <- 200
x <- 1:n # Must be equally spaced
phi <- 0.7
R <- phi^abs(outer(x, x, "-"))
sigma2 <- 1
set.seed(1330)
eta <- mvrnorm(1, rep(0, n), sigma2 * R)
cor(eta[1:(n - 1)], eta[2:n])
    ## [1] 0.7016662
    plot(x, eta, typ = "l")
abline(a = 0, b = 0)

    acf(eta)

    • Example correlation functions
      • Gaussian correlation function \[r_{ij}(\phi)=e^{-\frac{d_{ij}^{2}}{\phi}}\] where \(d_{ij}\) is the “distance” between locations i and j (note that \(d_{ij}=0\) for \(i=j\)) and \(r_{ij}(\phi)\) is the element in the \(i^{\textrm{th}}\) row and \(j^{\textrm{th}}\) column of \(\mathbf{R}(\phi)\).
    library(fields)
n <- 200
x <- 1:n
phi <- 40
d <- rdist(x)
R <- exp(-d^2/phi)
sigma2 <- 1

set.seed(4673)
eta <- mvrnorm(1, rep(0, n), sigma2 * R)
plot(x, eta, typ = "l")
abline(a = 0, b = 0)

    cor(eta[1:(n-1)], eta[2:n])
    ## [1] 0.9713954
    • Linear correlation function \[r_{ij}(\phi)=\begin{cases} 1-\frac{d_{ij}}{\phi} &\text{if}\ d_{ij}<0\\ 0 &\text{if}\ d_{ij}>0 \end{cases}\]
    library(fields)
n <- 200
x <- 1:n
phi <- 40
d <- rdist(x)
R <- ifelse(d<phi,1-d/phi,0)
sigma2 <- 1

set.seed(4803)
eta <- mvrnorm(1, rep(0, n), sigma2 * R)
plot(x, eta, typ = "l")
abline(a = 0, b = 0)

    cor(eta[1:(n-1)], eta[2:n])
    ## [1] 0.9775175

16.3 Bayesian Kriging

  • Read Chapter 17

  • Example data set

    url <- "https://www.dropbox.com/s/k2ajjvyxwvsbqgf/ISIT.txt?dl=1"
df <- read.table(url, header = TRUE)
df <- df[df$Station == 16,c(1,2,5,6)]

head(df)
    ##     Depth Sources Latitude Longitude
## 595   598   53.37  49.0322  -16.1493
## 596   631   51.32  49.0322  -16.1493
## 597   700   42.42  49.0322  -16.1493
## 598   666   38.32  49.0322  -16.1493
## 599  1317   37.29  49.0322  -16.1493
## 600  1352   36.27  49.0322  -16.1493
    plot(df$Depth, df$Sources, las = 1, ylim = c(0, 65), col = rgb(0, 0, 0, 0.25),
 pch = 19, xlab = "Depth (m)", ylab = "Sources")

  • Modify the linear model

    • Account for spatial autocorrelation
    • Enable prediction at locations
    • Write out on whiteboard
    • In R
    library(fields)
library(truncnorm)
library(mvtnorm)
library(coda)
library(MASS)

y <- df$Sources
X <- model.matrix(~I(Depth/Depth)-1,data=df)
Z <- diag(1,length(y))
n <- length(y)

# Distance matrix used to calculate correlation matrix
D <- rdist(df$Depth)

m.draws <- 10000  #Number of MCMC samples to draw
samples <- as.data.frame(matrix(,m.draws+1,dim(X)[2]+3))
colnames(samples) <- c(colnames(X),"sigma2.e","sigma2.alpha","phi")
samples[1,] <- c(40,10,150,250)  #Starting values for beta0, sigma2.e, sigma2.alpha, and phi
samples.alpha <- as.data.frame(matrix(,m.draws+1,length(y)))

# Hyperparameters for priors
sigma2.beta <- 100 
a <- 2.01 
b <- 1
c <- 2.01 
d <- 1
e <- 0
f <- 10^6

accept <- matrix(, m.draws + 1, 1) # Monitor acceptance rate
colnames(accept) <- c("phi")
phi.tune <- 10

set.seed(3909)
ptm <- proc.time()
for(k in 1:m.draws){
  beta <- t(samples[k,1:dim(X)[2]])
  sigma2.e <- samples[k,dim(X)[2]+1]
  sigma2.alpha <- samples[k,dim(X)[2]+2]
  phi <- samples[k,dim(X)[2]+3]
  C <- exp(-(D/phi)^2)
  CI <- solve(C)

  # Sample alpha
  H <- solve((1/sigma2.e)*t(Z)%*%Z + (1/sigma2.alpha)*CI)
  h <- (1/sigma2.e)*t(Z)%*%(y-X%*%beta)
  alpha <-  mvrnorm(1,H%*%h,H)

  # Sample beta
  H <- solve((1/sigma2.e)*t(X)%*%X + (1/sigma2.beta)*diag(1,dim(X)[2]))
  h <- (1/sigma2.e)*t(X)%*%(y-Z%*%alpha)
  beta <-  mvrnorm(1,H%*%h,H)

  # Sample sigma2.alpha
  sigma2.alpha <- 1/rgamma(1,a+n/2,b+t(alpha)%*%CI%*%(alpha)/2)

  # Sample sigma2.e
  sigma2.e <- 1/rgamma(1,c+n/2,d+t(y-X%*%beta-Z%*%alpha)%*%(y-X%*%beta-Z%*%alpha)/2)

  #Sample phi
  phi.star <- rnorm(1,phi,phi.tune)
  if(phi.star > e & phi.star < f){
C.star <-  exp(-(D/phi.star)^2)
mh1 <- dmvnorm(alpha,rep(0,n),sigma2.alpha*C.star,log=TRUE) + dunif(phi.star,e,f,log=TRUE)
mh2 <- dmvnorm(alpha,rep(0,n),sigma2.alpha*C,log=TRUE) + dunif(phi,e,f,log=TRUE)
R <- min(1,exp(mh1 - mh2))
if (R > runif(1)){phi <- phi.star}}

  # Save
  samples[k+1,] <- c(beta,sigma2.e,sigma2.alpha,phi)
  samples.alpha[k+1,] <- alpha
  accept[k+1,] <- ifelse(phi==phi.star,1,0)
  if(k%%1000==0) print(k)
}
proc.time() - ptm

mean(accept[-1,])

par(mfrow=c(1,1),mar=c(5,5,2,2))
plot(samples[,1],typ="l",xlab="k",ylab=expression(beta[0]))
plot(samples[,2],typ="l",xlab="k",ylab=expression(sigma[epsilon]^2))
plot(samples[,3],typ="l",xlab="k",ylab=expression(sigma[alpha]^2))
plot(samples[,4],typ="l",xlab="k",ylab=expression(phi))

#Note: there is a trace plots for each of the 51 alphas
plot(samples.alpha[,1],typ="l",xlab="k",ylab=expression(alpha))

burn.in <- 1000

#Effective sample size
effectiveSize(samples[-c(1:burn.in),])
effectiveSize(samples.alpha[-c(1:burn.in),])

# E() of beta, sigma2.e, sigma2.alpha, phi
colMeans(samples[-c(1:burn.in),])  

# 95% Equal-tailed CIs
apply(samples[-c(1:burn.in),],2,FUN=quantile,prob=c(0.025,0.975)) 


# obtain E(y) of at new depths
new.depth <- seq(min(df$Depth),max(df$Depth),by=10)
X.new <- model.matrix(~I(new.depth/new.depth)-1)
Z.new <- diag(1,length(new.depth))
D.cross <- rdist(new.depth,df$Depth)

samples.pred <- matrix(,length(burn.in:m.draws),length(new.depth))

for(k in burn.in:m.draws){
  alpha <- t(samples.alpha[k,])
  beta <- t(samples[k,1:dim(X)[2]])
  sigma2.e <- samples[k,dim(X)[2]+1]
  sigma2.alpha <- samples[k,dim(X)[2]+2]
  phi <- samples[k,dim(X)[2]+3]

  C <- exp(-(D/phi)^2) 
  C.cross <- exp(-(D.cross/phi)^2)

  alpha.new <- C.cross%*%solve(C)%*%alpha
  samples.pred[k-burn.in+1,] <- X.new%*%beta + Z.new%*%alpha.new
}


plot(df$Depth,df$Sources,las=1,ylim=c(0,60),col=rgb(0,0,0,0.25),pch=19,
 xlab="Depth (m)",ylab="Sources")
points(new.depth,colMeans(samples.pred),typ="l")


EV.beta <- mean(samples[-c(1:burn.in),1])  
EV.alpha <- colMeans(samples.alpha[-c(1:burn.in),])
EV.e <- y - X%*%EV.beta - Z%*%EV.alpha 

par(mfrow=c(2,1))  
plot(df$Depth, EV.e, las = 1, col = rgb(0, 0, 0, 0.25),pch = 19, xlab = "Depth (m)", ylab = expression("E("*epsilon*")"))
hist(EV.e,col="grey",main="",breaks=10,xlab= expression("E("*epsilon*")"))