Gp Nimble V2

#  devtools::install_github("cboettig/gpmanagement")

library("gpmanagement")
library("dplyr")
library("ggplot2")

  set.seed(1234)
  Tobs <- 40
  f <- function (x, h, p){
    sapply(x, function(x) {
        x <- pmax(0, x - h)
        x * exp(p[1] * (1 - x/p[2]) * (x - p[3])/p[2])
    })
  }
  p <- c(2, 8, 5)
  sigma_g <- 0.05
  z_g <- function() rlnorm(1, 0, sigma_g)
  x <- numeric(Tobs)
  x[1] <- 5.5
  for(t in 1:(Tobs-1))
    x[t+1] = z_g() * f(x[t], h=0, p=p)
  obs <- data.frame(x = c(0, 
                          pmax(rep(0,Tobs-1), x[1:(Tobs-1)])), 
                    y = c(0, 
                          x[2:Tobs]))
  xObs <- obs$x
  yObs <- obs$y
  xPred <- seq(0, 15, length=50)

  fit <- gp_setup(xObs, yObs, xPred)

[1] RW_block sampler: rho, sigGP, sigOE,  adaptive: TRUE,  adaptScaleOnly: FALSE,  adaptInterval: 200,  scale: 1,  propCov: identity

  Cmcmc <- fit$Cmcmc 
  Cpred <- fit$Cpred
  Cmodel <- fit$Cmodel
  
  system.time(Cmcmc$run(100000))

   user  system elapsed 
 25.548   0.024  25.734 

  samples <- as.matrix(Cmcmc$mvSamples)
  ## basic sanity check
  testthat::expect_identical(Cmodel$getNodeNames(topOnly = TRUE), colnames(samples))

predict from GP model using posterior MCMC samples

system.time(Cpred$run(samples))

   user  system elapsed 
 26.604   0.004  26.680 

extract predictions: E and C

  E <- Cpred$getE()
  C <- Cpred$getC()
  
obs <- data.frame(x = xObs, y = yObs)
pred <- data.frame(x = xPred, y = E, ymin = E - sqrt(diag(C)), ymax = E + sqrt(diag(C)))
ggplot2::ggplot(pred) + 
  geom_ribbon(aes(x = x,y = y, ymin = ymin, ymax = ymax), fill = "grey80") +
  geom_line(aes(x = x, y = y), size=1) + 
  geom_point(data = obs, aes(x,y)) +
  coord_cartesian(xlim = range(c(xObs, xPred)), ylim = range(c(yObs,E))) +
  theme_bw()