In cboettig/nonparametric-bayes: Nonparametric Bayes inference for ecological models

require(MASS)
require(ggplot2)
require(kernlab)
require(nonparametricbayes)
opts_knit$set(upload.fun = socialR::notebook.url)
opts_chunk$set(dev.args=list(bg="transparent"), comment=NA, tidy=FALSE, message=FALSE)
theme_update(plot.background = element_rect(fill = "transparent",colour = NA))

X <- seq(-5,5,len=50)
obs <- data.frame(x = c(-4, -3, -1,  0,  2, 3),
                  y = c(-2,  0,  1,  2, -1, -1))
l <- 1
sigma_n <- 0.1

#knit("beverton_holt_data.Rmd")
#true <- data.frame(x = X, y = sapply(X, f, 0, p))

Cholesky method

  SE <- function(Xi,Xj, l=1) exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2)
  cov <- function(X, Y) outer(X, Y, SE, l) 
  n <- length(obs$x)
  K <- cov(obs$x, obs$x)
  I <-  diag(1, n)
  L <- chol(K + sigma_n^2 * I)
  alpha <- forwardsolve(t(L), backsolve(L, obs$y)) # see http://stat.ethz.ch/R-manual/R-patched/library/base/html/chol2inv.html
  k_star <- cov(obs$x, X)
  Y <- t(k_star) %*% alpha
  v <- backsolve(L, k_star)
  Var <- cov(X,X) - t(v) %*% v
  loglik <- -.5 * t(obs$y) %*% alpha - sum(log(diag(L))) - n * log(2 * pi) / 2

Direct method

  d <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "direct")
  Ef <- d$Ef
  Cf <- d$Cf

Direct sequential method, avoids matrix inverse instability

s <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "sequential")
ef <- s$Ef
cf <- s$Cf

ggplot(data.frame(x=X, Ef=Ef, ef=ef))+ geom_point(aes(x,Ef), col='red') + geom_line(aes(x,ef))
````


### kernlab method

``` {r}
k <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "kernlab")
y_p <- k$Ef

Compare these results:

``` {r} require(reshape2) df <- data.frame(x = X, direct = Ef, Cholesky = Y, kernlab = y_p, sequential = ef) df <- melt(df, id = "x") ggplot(df)+ geom_jitter(aes(x, value, color = variable)) + geom_point(data = obs, aes(x,y)) ````

cboettig/nonparametric-bayes documentation built on May 13, 2019, 2:09 p.m.