object$meanResponse <- added
}
return(object)
}
myem.exp1 <- EMULATElm.GPEm(exp1)
AddBest
eval(parse(text=paste("lm(",exp1$Response,"~1,data=",exp1$dString,")",sep="")))
exp1$Response
exp1$sString
exp1$dString
exp1 <- GPEm(tData = tData, Response = names(tData)[4], tcands = cands,
tcanfacs = NULL, dString = "tData", TryFouriers = TRUE, maxOrder = 2,
maxdf = 2)
exp1
exp1$Response
GPEm <- function(tData = NULL, Response = NULL, XX = NULL, tcands, tcanfacs = NULL,
dString = "tData",
meanResponse = NULL, nugget = NULL, maxdf = NULL, TryFouriers = FALSE,
maxOrder = NULL, additionalVariables = NULL, params = NULL,
Zp.mean = NULL, Zp.var = NULL, ZZ.mean = NULL, ZZ.var = NULL) {
out <- list(X = tData[tcands], n = nrow(tData[tcands]), d = ncol(tData[tcands])-1,
Z = tData$Response, Response = Response, XX = XX, nn = nrow(XX), meanResponse = meanResponse,
tcands = tcands, tcanfacs = tcanfacs, nugget = nugget, maxdf = maxdf,
TryFouriers = TryFouriers, maxOrder = maxOrder, dString = dString,
additionalVariables = additionalVariables, params = params,
Zp.mean = Zp.mean, Zp.var = Zp.var, ZZ.mean = ZZ.mean, ZZ.var = ZZ.var)
class(out) <- "GPEm"
return(out)
}
exp1 <- GPEm(tData = tData, Response = names(tData)[4], tcands = cands,
tcanfacs = NULL, dString = "tData", TryFouriers = TRUE, maxOrder = 2,
maxdf = 2)
exp1$Response
myem.exp1 <- EMULATElm.GPEm(exp1)
tData
library(sensitivity)
?tell
?sobolroalhs
x <- sobolroalhs(model = NULL, factors = 3, N = 1000, order =1, nboot=0)
x
x$X <- 2*randomLHS(100,3) - 1
x$X
x <- sobolroalhs(model = NULL, factors = 3, N = 1000, order =1, nboot=0)
x$X <- 2*x$X - 1
library(ncdf4)
rm(list=ls())
#source('~/Dropbox/CCCMA/2016/DannyDevelopment.R')
twd <- getwd()
#setwd("~/Dropbox/French Tuning/POUR_DANIEL/SCM")
setwd("/Users/victoriavolodina/Desktop/French_meeting/POUR_DANIEL/SCM")
nc000 <- nc_open(file="Out_klevel_REF629GCM.tuning.01.nc")
#nc000
levf <- ncvar_get(nc000,"levf")
time <- ncvar_get(nc000,"time")
pTemp <- ncvar_get(nc000,"theta")
#pTemp
th_bl <- ncvar_get(nc000,"th_bl")
days <- time/(3600*24)
tfile <- "Out_klevel_REF629GCM.tuning"
alts <- paste("0",as.character(1:9),sep="")
alts <- c(alts,"10")
tfiles <- paste(tfile,alts,"nc",sep=".")
th_bls <- matrix(NA,nrow=length(tfiles),ncol=length(th_bl))
for(k in 1:length(tfiles)){
ncdummy <- nc_open(file=tfiles[k])
thDummy <- ncvar_get(ncdummy,"th_bl")
th_bls[k,] <- thDummy
}
p
parameters <- read.table("simulations.txt",colClasses = "character")
RKDX = parameters$V1
RKDX = unlist(lapply(strsplit(RKDX,"="), function(e) e[2]))
TENTRX = parameters$V2
TENTRX = unlist(lapply(strsplit(TENTRX,"="), function(e) e[2]))
ALFX = parameters$V3
ALFX = unlist(lapply(strsplit(ALFX,"="), function(e) e[2]))
tParams <- data.frame(RKDX=as.numeric(RKDX),TENTRX=as.numeric(TENTRX),ALFX=as.numeric(ALFX))
tmins <- apply(tParams,2,min)
tmaxs <- apply(tParams,2,max)
tRanges <- tmaxs-tmins
scaledParams <- sapply(1:3, function(i) 2*((tParams[,i] - tmins[i])/tRanges[i]) - 1)
tData <- data.frame(cbind(scaledParams,th_bls[,13]))
names(tData) <- c(names(tParams),"thbls13")
Noise <- rnorm(10,mean(tData[,4]),sd(tData[,4]))
tData <- cbind(tData,Noise)
cands <- names(tData)[c(1:3,5)]
setwd("/Users/victoriavolodina/Desktop/French_meeting")
source("~/Desktop/French_meeting/StanEmulateCode_VVMyVersion.R" )
myem.lm <- EMULATE.lm(Response=names(tData)[4], tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 2)
myem.gp = EMULATE.gpstan(meanResponse=myem.lm, nugget=0.0001, tData=tData, additionalVariables=names(tData)[1:3],CreatePredictObject = TRUE)
EMULATOR.gpstan <- function(x, Emulator, GP=TRUE){
if(!GP){# || (Emulator$nugget==1)){
lin <- Emulator$linModel
newDat <- as.data.frame(x)
prediction <- try(custom.predict(lin,newdata=newDat,Emulator$pre.Lists),silent=TRUE)
if(inherits(prediction,"try-error")){
stop("Ensure the data frame has the correct named columns")
}
return(list(Expectation=as.vector(prediction[,1]), Variance = as.vector(prediction[,2])^2))
}
else{
Xpred <- as.data.frame(x)
tt = terms(Emulator$linModel)
Terms = delete.response(tt)
mm = model.frame(Terms, Xpred, xlev = Emulator$linModel$xlevels)
Hpred = model.matrix(Terms, mm, contrasts.arg = Emulator$linModel$contrasts)
if(is.null(Emulator$StanPredict))
fit.y2 <- stan(file = tprednewfile_loc,
data = list(X1 = Emulator$Design, X2 = Xpred, y1= Emulator$tF,
H1 = Emulator$H, H2 = Hpred, N1 = dim(Emulator$Design)[1], N2 = dim(Xpred)[1],
Np = dim(Emulator$H)[2], p = dim(Emulator$Design)[2], M = dim(Emulator$ParameterSamples$beta)[1],
sigma_sq = Emulator$ParameterSamples$sigma_sq, beta = Emulator$ParameterSamples$beta,
delta = Emulator$ParameterSamples$delta),
iter = 1, warmup = 0, chains = 1, cores = 1, pars = c("tmeans","tsds"), include=TRUE,
algorithm = c('Fixed_param'))
else
fit.y2 <- stan(fit=Emulator$StanPredict,
data = list(X1 = Emulator$Design, X2 = Xpred, y1= Emulator$tF,
H1 = Emulator$H, H2 = Hpred, N1 = dim(Emulator$Design)[1], N2 = dim(Xpred)[1],
Np = dim(Emulator$H)[2], p = dim(Emulator$Design)[2], M = dim(Emulator$ParameterSamples$beta)[1],
sigma_sq = Emulator$ParameterSamples$sigma_sq, beta = Emulator$ParameterSamples$beta,
delta = Emulator$ParameterSamples$delta),
iter = 1, warmup = 0, chains = 1, cores = 1, pars = c("tmeans","tsds"), include=TRUE,
algorithm = c('Fixed_param'))
predict.y2 <- extract(fit.y2, pars = c('tmeans','tsds'))
tExpectation <- predict.y2$tmeans[1,]
StandardDev <- predict.y2$tsds[1,]
predict.y2 <- summary(fit.y2, pars = c('tmeans', 'tsds'), use_cache = FALSE)$summary[, 1]
# NN <- length(predict.y2)
#  tExpectation <- as.numeric(predict.y2[1:NN/2])
#  StandardDev <- as.numeric(predict.y2[(NN/2+1):NN])
return(list(Expectation=tExpectation,Variance=StandardDev^2))
}
}
EMULATE.lm <- function(Response, tData, dString="tData",maxdf=NULL,tcands=cands,tcanfacs=canfacs,TryFouriers=FALSE,maxOrder=NULL){
if(is.null(maxdf)){
maxdf <- ceiling(length(tData[,1])/10)
}
startLM <- eval(parse(text=paste("lm(",Response,"~1,data=",dString,")",sep="")))
if(TryFouriers)
msl <- list(linModel=startLM,Names=NULL,mainEffects=NULL,Interactions=NULL,Factors=NULL,FactorInteractions=NULL,ThreeWayInters=NULL,DataString=dString,ResponseString=Response,tData=tData,BestFourier=TRUE,maxOrder=maxOrder)
else
msl <- list(linModel=startLM,Names=NULL,mainEffects=NULL,Interactions=NULL,Factors=NULL,FactorInteractions=NULL,ThreeWayInters=NULL,DataString=dString,ResponseString=Response,tData=tData,BestFourier=FALSE)
added <- AddBest(tcands,tcanfacs,msl)
for(i in 1:30){
added <- AddBest(tcands,tcanfacs,added)
if(!is.null(added$Break))
break
}
print(summary(added$linModel))
if(length(tData[,1])-1-added$linModel$df > maxdf){
trm <- removeNterms(N=500,linModel=added$linModel,dataString=added$DataString,responseString=added$ResponseString,Tolerance=sum(sort(anova(added$linModel)$"Sum Sq"))*(1e-4),Names=added$Names,mainEffects=added$mainEffects,Interactions=added$Interactions,Factors=added$Factors,FactorInteractions=added$FactorInteractions,ThreeWayInters=added$ThreeWayInters,tData=added$tData,Fouriers = added$Fouriers)
print(summary(trm$linModel))
trm2 <- removeNterms(N=max(c(0,length(tData[,1])-maxdf-1-trm$linModel$df)),linModel=trm$linModel,dataString=added$DataString,responseString=added$ResponseString,Tolerance=sum(sort(anova(added$linModel)$"Sum Sq"))*(5e-1),Names=trm$Names,mainEffects=trm$mainEffects,Interactions=trm$Interactions,Factors=trm$Factors,FactorInteractions=trm$FactorInteractions,ThreeWayInters=trm$ThreeWayInters,tData=added$tData,Fouriers=trm$Fouriers)
print(summary(trm2$linModel))
trm2$pre.Lists <- get.predict.mats(trm2$linModel)
trm2$DataString <- dString
trm2$ResponseString <- Response
return(trm2)
}
else{
return(added)
}
}
EMULATE.gpstan <- function(meanResponse, nugget, tData, additionalVariables=NULL, CreatePredictObject=FALSE, ...){
#additionalVariables is a vector of parameter names that we want to fit a GP to that didn't make it into the model
#the default model, with additionalVariables=NULL, is to only treat those terms that make it into the linear model as active
#tData is a data frame containing the inputs and outputs and must be the same as that used to create meanResponse
Design <- as.matrix(tData[,which((names(tData)%in%meanResponse$Names) | (names(tData)%in%meanResponse$Factors) | (names(tData)%in%additionalVariables) | (names(tData) %in% names(meanResponse$Fouriers)))])
tF <- tData[,which(names(tData)==meanResponse$ResponseString)]
H <- model.matrix(meanResponse$linModel)
consEm <- EMULATE.lm(Response=meanResponse$ResponseString, tData=tData, tcands="Noise",tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 0)
sd2 <- (summary(consEm$linModel)$sig - summary(meanResponse$linModel)$sig)^2
sigsq <- summary(meanResponse$linModel)$sigma^2
sigsqvar <- (sd2)^2
StanEmulator <- stan(file=tfile_loc,
data=list(X1 = Design, y1= tF, H1 = H,  N1 = dim(Design)[1], Np = dim(H)[2], p = dim(Design)[2],
SigSq = sigsq, SigSqV = sigsqvar, nugget=nugget),
iter = 4000,warmup=3000, chains = 2, cores = 4, ...)
ParameterSamples <- extract(StanEmulator, pars = c('sigma_sq', 'delta_par', 'beta'))
#To avoid recompiling in EMULATOR.gpstan, create stan prediction object
if(CreatePredictObject){
Xpred <- 2*randomLHS(10,length(Design[1,])) - 1
Xp <- as.data.frame(Xpred)
names(Xp) <- names(tData)[1:length(Design[1,])]
tt = terms(meanResponse$linModel)
Terms = delete.response(tt)
mm = model.frame(Terms, Xp, xlev = meanResponse$linModel$xlevels)
Hpred = model.matrix(Terms, mm, contrasts.arg = meanResponse$linModel$contrasts)
StanPredict <- stan(file = tprednewfile_loc,
data = list(X1 = Design, X2 = Xp, y1= tF,
H1 = H, H2 = Hpred, N1 = dim(Design)[1], N2 = dim(Xp)[1],
Np = dim(H)[2], p = dim(Design)[2], M = dim(ParameterSamples$beta)[1],
sigma_sq = ParameterSamples$sigma_sq, beta = ParameterSamples$beta,
delta = ParameterSamples$delta),
iter = 1, warmup = 0, chains = 1, cores = 1, pars = c("tmeans","tsds"), include=TRUE,
algorithm = c('Fixed_param'))
}
else
StanPredict <- NULL
gp.list <- list(Design=Design, tF=tF, H=H, ParameterSamples=ParameterSamples, StanPredict=StanPredict)
return(c(meanResponse,gp.list))
}
ValidPlot <- function(fit, X, y, interval = c(), axis = 1, heading = " ", xrange = c(-1, 1), ParamNames) {
# Plots the emulator predictions together with error bars and true values.
#
# Args:
#  fit: A matrix of emulator predictions. 1st column corresponds to
#             posterior mean values. 2nd column to posterior mean minus 2
#             * standard deviation and 3rd column to posterior mean plus 2
#             * standard deviation.
#  X: A matrix of inputs.
#  y: Simulator evaluated at X.
#  interval: A vector to define the y limit of the plot (optional).
#  axis: Input against which the predictions are plotted.
#  heading: A heading for the plot.
#  xrange:  A vector to define the x limit of the plot.
#  ParamNames: A vector of names of the parameters in the data
#
# Returns:
#  A plot of emulator predictions together with error bars
#  (plus and minus 2* standard deviations) and true values.
outside <- c()
inside <- c()
maxval <- max(fit)
minval <- min(fit)
for(i in 1:length(y)) {
if(fit[i, 2] <= y[i] & y[i] <= fit[i, 3]) {
inside <- c(inside,i)
}
}
outside <- c(1:length(y))[-inside]
if(is.null(interval)) {
plot(X[ , axis], fit[, 1], pch=20, ylab='Y',
xlab=ParamNames[axis], ylim=c(minval, maxval),
xlim=xrange, main=heading,cex.main=0.8)
} else {
plot(X[ , axis], fit[, 1], pch=20, ylab='Y',
xlab=ParamNames[axis], ylim=interval,
xlim=xrange, main=heading,cex.main=0.8)
}
arrows(X[, axis], fit[, 2], X[, axis], fit[, 3],
length=0.05, angle=90, code=3, col='black')
points(X[inside, axis], y[inside], pch=20, col='green')
points(X[outside, axis], y[outside], pch=20, col='red')
}
myem.gp = EMULATE.gpstan(meanResponse=myem.lm, nugget=0.0001, tData=tData, additionalVariables=names(tData)[1:3],CreatePredictObject = TRUE)
myem.lm
stanLOOplot(StanEmulator = myem.gp, ParamNames=names(tData)[1:3])
myem.gp = EMULATE.gpstan(meanResponse=myem.lm, nugget=0.0001, tData=tData, additionalVariables=names(tData)[1:3],CreatePredictObject = TRUE)
tData
library(lhs)
library(rstan)
source('/Users/victoriavolodina/Desktop/Banerjee Method Paper/Academic Example/Functions.R')
source("~/Desktop/NewStan/AutoLMcodeMyVersion.R")
source("~/Desktop/NewStan/CustomPredictMyVersion.R")
source("~/Desktop/NewStan/impLayoutplotMyVersion.R")
source("~/Desktop/NewStan/StanEmulateCode_VVMyVersion.R")
ToyFun <- function(X) {
return(sin(1/((X[, 1]*0.7+0.3)*(X[, 2]*0.7+0.3))))
}
x1<-c(0,.02,.075,.08,.14,.15,.155,.156,.18,.22,.29,.32,.36,
.37,.42,.5,.57,.63,.72,.785,.8,.84,.925,1)
x2<-c(.29,.02,.12,.58,.38,.87,.01,.12,.22,.08,.34,.185,.64,
.02,.93,.15,.42,.71,1,0,.21,.5,.785,.21)
X.D1 <- cbind(x1, x2)
Y.D1 <- ToyFun(X.D1)
xx <- seq(0, 1, length.out=30)
grid <- expand.grid(x=xx, y=xx)
X.valid.grid <- as.matrix(grid, nrow=900, ncol=2)
WavyValues <- matrix(ToyFun(grid), nrow=30, ncol=30)
quartz()
par(mfrow=c(1, 1))
res <- persp(xx, xx, WavyValues,
main='Wavy', xlab='X1', ylab='X2', zlab='Y',
theta=35, phi=30, col='lightblue')
mypoints <- trans3d(X.D1[, 1], X.D1[, 2], Y.D1, pmat = res)
points(mypoints, pch=20, col='red')
X.S <- matrix(c(0.8, 0.75), ncol=2)
Z <- ToyFun(X.S) + rnorm(1, 0, 0.05) + rnorm(1, 0, 0.1)
mypoints <- trans3d(X.S[, 1], X.S[, 2], Z, pmat = res)
points(mypoints, pch=8, col='yellow')
H1 <- cbind(rep(1, dim(X.D1)[1]), X.D1) # construct regression matrix
fit1.exp <- stan(file = '/Users/victoriavolodina/Desktop/NewStan/MySpeed1MyVersion.stan',
data = list(X1 = X.D1, y1 = Y.D1, N1 = 24, Np = 3, p = 2,
H1 = H1, nugget = 0.0001), iter = 2000, chains = 2, cores = 4)
summary(fit.exp1)
names(fit1.exp)
extract(fit1.param$lp_)
extract(fit1.exp$lp_)
extract(fit1.exp, pars = c('lp_'))
extract(fit1.exp, pars = c('lp__'))
hist(extract(fit1.exp, pars = c('lp__')))
lp <- extract(fit1.exp, pars = c('lp__'))
lp
dim(lp)
length(lp)
hist(lp[[1]])
fit1.param <- extract(fit1.exp, pars = c('sigma_sq', 'delta_par', 'beta'))
plot(fit1.param$sigma_sq, lp[[1]])
plot(fit1.param$sigma_sq, exp(lp[[1]]))
quartz()
plot(fit1.param$sigma_sq, lp[[1]])
library(lhs)
library(rstan)
source('/Users/victoriavolodina/Desktop/Banerjee Method Paper/Academic Example/Functions.R')
source("~/Desktop/NewStan/AutoLMcodeMyVersion.R")
source("~/Desktop/NewStan/CustomPredictMyVersion.R")
source("~/Desktop/NewStan/impLayoutplotMyVersion.R")
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
ToyFun <- function(X) {
return(sin(1/((X[, 1]*0.7+0.3)*(X[, 2]*0.7+0.3))))
}
x1<-c(0,.02,.075,.08,.14,.15,.155,.156,.18,.22,.29,.32,.36,
.37,.42,.5,.57,.63,.72,.785,.8,.84,.925,1)
x2<-c(.29,.02,.12,.58,.38,.87,.01,.12,.22,.08,.34,.185,.64,
.02,.93,.15,.42,.71,1,0,.21,.5,.785,.21)
Y.D1 <- ToyFun(cbind(x1, x2))
X <- 2*cbind(x1, x2)-1
tData.toy <- as.data.frame(cbind(X, Y.D1))
names(tData.toy) <- c('x1', 'x2', 'Y')
toy.lm <- EMULATE.lm(Response=names(tData.toy)[3], tData=tData.toy,
tcands=names(tData.toy)[c(1, 2)],tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 2)
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = TRUE)
names(toy.gp)
names(toy.gp$ParameterSamples)
length(toy.gp$ParameterSamples)
toy.gp$ParameterSamples[[1]]
toy.gp$ParameterSamples$sigma_sq
?list
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = TRUE)
names(toy.gp)
toy.gp$ParameterSamples
toy.gp$ParameterSamples$sigma_sq
stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
toy.gp$ParameterSamples$beta
dim(toy.gp$ParameterSamples$beta)
matrix(toy.gp$ParameterSamples$beta, ncol=1)
matrix(toy.gp$ParameterSamples$beta, nrow=1)
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = TRUE)
names(toy.gp)
toy.gp$ParameterSamples
stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
length(toy.gp$ParameterSamples$sigma_sq)
dim(toy.gp$ParameterSamples$beta)[1]
tpredfile_loc
stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
toy.LOO <- stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
names(toy.LOO)
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = TRUE)
names(toy.gp)
toy.gp$ParameterSamples
toy.LOO <- stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
names(toy.gp)
dim(toy.gp$ParameterSamples$sigma_sq)
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = TRUE)
names(toy.gp)
toy.gp$ParameterSamples
toy.LOO <- stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
toy.gp$ParameterSamples$sigma_sq
toy.gp$ParameterSamples$sigma_sq[1, ]
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = TRUE)
library(lhs)
library(rstan)
source('/Users/victoriavolodina/Desktop/Banerjee Method Paper/Academic Example/Functions.R')
#source('/Users/victoriavolodina/Desktop/History Matching/HMFun.R')
source("~/Desktop/NewStan/AutoLMcodeMyVersion.R")
source("~/Desktop/NewStan/CustomPredictMyVersion.R")
source("~/Desktop/NewStan/impLayoutplotMyVersion.R")
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
ToyFun <- function(X) {
return(sin(1/((X[, 1]*0.7+0.3)*(X[, 2]*0.7+0.3))))
}
x1<-c(0,.02,.075,.08,.14,.15,.155,.156,.18,.22,.29,.32,.36,
.37,.42,.5,.57,.63,.72,.785,.8,.84,.925,1)
x2<-c(.29,.02,.12,.58,.38,.87,.01,.12,.22,.08,.34,.185,.64,
.02,.93,.15,.42,.71,1,0,.21,.5,.785,.21)
Y.D1 <- ToyFun(cbind(x1, x2))
X <- 2*cbind(x1, x2)-1
tData.toy <- as.data.frame(cbind(X, Y.D1))
names(tData.toy) <- c('x1', 'x2', 'Y')
toy.lm <- EMULATE.lm(Response=names(tData.toy)[3], tData=tData.toy,
tcands=names(tData.toy)[c(1, 2)],tcanfacs=NULL,TryFouriers=TRUE,maxOrder=2,maxdf = 2)
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = FALSE)
AA <- matrix(c(2, 3, 1, 15, 7, 4, 4, 19, 9), nrow=3, ncol=3)
AA
Delta <- diag(c(0.5, 4, 1.5))
Delta
AA*Delta
Delta*AA
Delta*t(AA)
Xpred <- randomLHS(10, 2)
Xp <- as.data.frame(2*Xpred-1)
names(Xp) <- names(tData.toy)[1:2]
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
names(NewSamples)
NewSamples$Expectation
Xpred <- randomLHS(100, 2)
Xp <- as.data.frame(2*Xpred-1)
names(Xp) <- names(tData.toy)[1:2]
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples1 <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples1 <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
if(!file.exists("~/.R/Makevars")) file.show("~/.R/Makevars")
names(toy.gp)
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
tfile_loc
system.time(toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = FALSE))
system.time(toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = FALSE))
tfile_loc = "~/Desktop/NewStan/MySpeed1mod.stan"
system.time(toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = FALSE))
toy.LOO <- stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
Xpred <- randomLHS(100, 2)
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
tprednewfile_loc
Xpred <- randomLHS(100, 2)
Xp <- as.data.frame(2*Xpred-1)
names(Xp) <- names(tData.toy)[1:2]
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
names(NewSamples)
fitt <- cbind(NewSamples$Expectation, NewSamples$Expectation - 2*NewSamples$Expectation,
NewSamples$Expectation + 2*NewSamples$Expectation)
fitt
valid.NewSamples <- ToyFun(Xp)
quartz()
ValidPlot(fitt, Xp, valid.NewSamples, interval = c(), axis = 1, heading = " ", xrange = c(-1, 1), ParamNames=c('x1', 'x2'))
valid.NewSamples <- ToyFun(Xpred)
quartz()
ValidPlot(fitt, Xp, valid.NewSamples, interval = c(), axis = 1, heading = " ", xrange = c(-1, 1), ParamNames=c('x1', 'x2'))
fitt <- cbind(NewSamples$Expectation, NewSamples$Expectation - 2*sqrt(NewSamples$Variance),
NewSamples$Expectation + 2*sqrt(NewSamples$Variance))
quartz()
ValidPlot(fitt, Xp, valid.NewSamples, interval = c(), axis = 1, heading = " ", xrange = c(-1, 1), ParamNames=c('x1', 'x2'))
NewSamples1 <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
NewSamples1 <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
fitt <- cbind(NewSamples1$Expectation, NewSamples1$Expectation - 2*sqrt(NewSamples1$Variance),
NewSamples1$Expectation + 2*sqrt(NewSamples1$Variance))
quartz()
ValidPlot(fitt, Xp, valid.NewSamples, interval = c(), axis = 1, heading = " ", xrange = c(-1, 1), ParamNames=c('x1', 'x2'))
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
fitt <- cbind(NewSamples$Expectation, NewSamples$Expectation - 2*sqrt(NewSamples$Variance),
NewSamples$Expectation + 2*sqrt(NewSamples$Variance))
quartz()
ValidPlot(fitt, Xp, valid.NewSamples, interval = c(), axis = 1, heading = " ", xrange = c(-1, 1), ParamNames=c('x1', 'x2'))
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
rfile_loc
tfile_loc
tpredfile_loc
tprednewfile_loc
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = FALSE)
toy.LOO <- stanLOOplot(StanEmulator = toy.gp, ParamNames=names(tData.toy)[1:2])
NewSamples <- EMULATOR.gpstan(Xp, Emulator=toy.gp)
fitt <- cbind(NewSamples$Expectation, NewSamples$Expectation - 2*sqrt(NewSamples$Variance),
NewSamples$Expectation + 2*sqrt(NewSamples$Variance))
quartz()
ValidPlot(fitt, Xp, valid.NewSamples, interval = c(), axis = 1, heading = " ", xrange = c(-1, 1), ParamNames=c('x1', 'x2'))
source("~/Desktop/NewStan/StanEmulateCodeMyVersion.R")
tfile_loc
toy.gp <- EMULATE.gpstan(meanResponse=toy.lm, nugget=0.0001, tData=tData.toy, additionalVariables=NULL,CreatePredictObject = FALSE,
FastVersion = FALSE)
if (!file.exists("~/.R/Makevars")) {
cat('CXX=g++ -arch x86_64 -ftemplate-depth-256 -stdlib=libstdc++\n
CXXFLAGS="-mtune=native  -O3 -Wall -pedantic -Wconversion"\n',
file="~/.R/Makevars");
} else {
file.show("~/.R/Makevars");
}
library(ncdf4)
source('~/Desktop/NewStan/StanEmulateCodeRMyVersion.R')
setwd("~/Documents/ExeterUQ")
load("/Users/victoriavolodina/Documents/GitFrench/cnrmWave1.RData")
scaledParams <- load("/Users/victoriavolodina/Documents/GitFrench/cnrmMasterWave1_US.RData")
scaledParams <- cnrmMasterWave1_US
FullData <- tData
tData <- tData[1:120, ]
cands <- c(names(scaledParams)[-c(1, 11:14)], "Noise")
par(mfrow = c(3, 3), mar = c(4, 4, 2, 1))
for(i in 1:(length(cands)-1)) plot(tData[, cands[i]], tData[, "SCMdata"], pch = 20, xlab = cands[i],
ylab = 'Pot Temp Average')
myem.lm.cnrm <- EMULATE.lm(Response="SCMdata", tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=4,maxdf = 30)
source('BuildEmulator/BuildEmulator.R')
myem.lm.cnrm <- EMULATE.lm(Response="SCMdata", tData=tData, tcands=cands,tcanfacs=NULL,TryFouriers=TRUE,maxOrder=4,maxdf = 30)
myem.lm.cnrm$linModel <- lm(SCMdata ~ ALFX + TENTRX+TENTR + RKDX+VVX +VVN+GAMAP1+REFLCAPE, data = tData)
myem.gp.cnrm.nugget <- EMULATE.gpstan(meanResponse=myem.lm.cnrm, nugget = 0.0001, nuggetPrior = TRUE, tData=tData, additionalVariables=NULL,CreatePredictObject = FALSE,FastVersion = TRUE)