########## R script: wolfcampBivPenSplLOOana.Rs ########## # For performing MCMC-based fitting of a bivariate penalized # spline fit to the Wolfcamp aquifer data. # Last changed: 28 MAY 2019 # Set tasks: doMCMC <- TRUE ; compileCode <- TRUE doLOOcalc <- TRUE options(warn=-1) # Load required packages: library(rstan) ; library(loo) # Load required function: source("Ztps.r") # Set MCMC parameters: nWarm <- 1000 # Length of warmup. nKept <- 1000 # Size of the kept sample. nThin <- 1 # Thinning factor. # Read in and extract data: wolfcamp <- read.table("wolfcamp.txt",header=TRUE) yOrig <- wolfcamp$piezometricHead x1Orig <- wolfcamp$xCoordinate x2Orig <- wolfcamp$yCoordinate # Standardised data for MCMC: y <- 2*(yOrig - min(yOrig))/(max(yOrig) - min(yOrig)) x1 <- 2*(x1Orig - min(x1Orig))/(max(x1Orig) - min(x1Orig)) x2 <- 2*(x2Orig - min(x2Orig))/(max(x2Orig) - min(x2Orig)) # Set up matrices for bivariate penalized spline model: wolfcampKnots <- read.table("wolfcampKnots.txt") X <- cbind(1,x1,x2) Z <- Ztps(cbind(x1,x2),wolfcampKnots) if (doMCMC) { # Specify model in Stan: bivPenSplModel <- 'data { int n; int ncX; int ncZ; vector[n] y; matrix[n,ncX] X; matrix[n,ncZ] Z; real sigmabeta; real ssigma; } parameters { vector[ncX] beta; real sigmaeps; vector[ncZ] u; real sigmau; } model { y ~ normal(X*beta + Z*u,sigmaeps); u ~ normal(0,sigmau); beta ~ normal(0,sigmabeta); sigmau ~ cauchy(0,ssigma); sigmaeps ~ cauchy(0,ssigma); } generated quantities { vector[n] log_lik; for (i in 1:n) { log_lik[i] = normal_lpdf(y[i]|X[i]*beta + Z[i]*u,sigmaeps); } }' # Store data in a list in format required by Stan: allData <- list(n=length(y),ncX=ncol(X),ncZ=ncol(Z),y=y,X=X, Z=Z,sigmabeta=1e5,ssigma=1e5) # Compile code for model if required: if (compileCode) stanCompilObj <- stan(model_code=bivPenSplModel,data=allData, iter=1,chains=1) # Obtain MCMC samples for each parameter using Stan: stanObj <- stan(model_code=bivPenSplModel,data=allData, warmup=nWarm,iter=(nWarm+nKept), chains=1,thin=nThin,refresh=250,fit=stanCompilObj) } if (doLOOcalc) { library(loo) LOOobj <- loo(extract_log_lik(stanObj)) cat("\n\n\n The leave-one-out information criterion for the bivariate spline model is:",LOOobj$estimates[3,1]) cat("\n\n\n") } ########## End of wolfcampBivPenSplLOOana ##########