Nothing
### R code from vignette source 'calex.Rnw'
###################################################
### code chunk number 1: LoadLibrary
###################################################
###################################################
### code chunk number 2: calex.Rnw:163-164
###################################################
library("calibrator")
###################################################
### code chunk number 3: size_of_datasets
###################################################
do_from_scratch = FALSE
n1 <- 20
n2 <- 21
###################################################
### code chunk number 4: make_D1
###################################################
D1.int <- latin.hypercube(n1,4)
rownames(D1.int) <- paste("coderun",1:nrow(D1.int),sep=".")
colnames(D1.int) <- c("x", "y", "A", "B")
head(D1.int)
###################################################
### code chunk number 5: make_D2
###################################################
D2.int <- latin.hypercube(n2,2)
rownames(D2.int) <- paste("obs",1:nrow(D2.int),sep=".")
colnames(D2.int) <- c("x", "y")
head(D2.int)
###################################################
### code chunk number 6: make_extractor.int
###################################################
extractor.int <- function(D1){
return(list(x.star = D1[, 1:2, drop = FALSE], t.vec = D1[, 3:4, drop = FALSE]))
}
###################################################
### code chunk number 7: make_h1.int
###################################################
h1.int <- function(xin){
out <- c(1,xin[1],xin[3],xin[1]*xin[3], xin[2]^2*xin[4])
names(out) <- c("const" , "x", "A", "Ax" , "By.sq")
return(out)
}
###################################################
### code chunk number 8: make_H1.int
###################################################
H1.int <- function(D1){
if (is.vector(D1)) {
D1 <- t(D1)
}
out <- t(apply(D1, 1, h1.int))
colnames(out) <- c("const" , "x", "A", "Ax" , "By.sq")
return(out)
}
###################################################
### code chunk number 9: make_h2.int
###################################################
h2.int <- function(x){
out <- c(x[1],x[1]*x[2])
names(out) <- c("h2.x" ,"h2.xy")
return(out)
}
H2.int <-
function (D2)
{
if (is.vector(D2)) {
D2 <- t(D2)
}
out <- t(apply(D2, 1, h2.int))
colnames(out) <- names(h2.int(D2[1, , drop = TRUE]))
return(out)
}
###################################################
### code chunk number 10: make_pdm.maker.psi1
###################################################
pdm.maker.psi1 <- function(psi1) {
jj.omega_x <- diag(psi1[1:2])
rownames(jj.omega_x) <- names(psi1[1:2])
colnames(jj.omega_x) <- names(psi1[1:2])
jj.omega_t <- diag(psi1[3:4],ncol=2)
rownames(jj.omega_t) <- names(psi1[3:4])
colnames(jj.omega_t) <- names(psi1[3:4])
sigma1squared <- psi1[5]
return(list(omega_x = jj.omega_x, omega_t = jj.omega_t,
sigma1squared = sigma1squared))
}
pdm.maker.psi2 <- function(psi2) {
jj.omegastar_x <- diag(psi2[1:2],ncol=2)
sigma2squared <- psi2[3]
return(list(omegastar_x = jj.omegastar_x, sigma2squared = sigma2squared))
}
###################################################
### code chunk number 11: phi.fun.int
###################################################
phi.fun.int <-
function (rho, lambda, psi1, psi1.apriori, psi2, psi2.apriori,
theta.apriori)
{
# CHANGES START
pdm.maker.psi1 <- function(psi1) {
jj.omega_x <- diag(psi1[1:2])
rownames(jj.omega_x) <- names(psi1[1:2])
colnames(jj.omega_x) <- names(psi1[1:2])
jj.omega_t <- diag(psi1[3:4],ncol=2)
rownames(jj.omega_t) <- names(psi1[3:4])
colnames(jj.omega_t) <- names(psi1[3:4])
sigma1squared <- psi1[5]
return(list(omega_x = jj.omega_x, omega_t = jj.omega_t,
sigma1squared = sigma1squared))
}
pdm.maker.psi2 <- function(psi2) {
jj.omegastar_x <- diag(psi2[1:2],ncol=2)
sigma2squared <- psi2[3]
return(list(omegastar_x = jj.omegastar_x, sigma2squared = sigma2squared))
}
# CHANGES END: remainder of function unaltered
jj.mean <- theta.apriori$mean
jj.V_theta <- theta.apriori$sigma
jj.discard.psi1 <- pdm.maker.psi1(psi1)
jj.omega_t <- jj.discard.psi1$omega_t
jj.omega_x <- jj.discard.psi1$omega_x
jj.sigma1squared <- jj.discard.psi1$sigma1squared
jj.discard.psi2 <- pdm.maker.psi2(psi2)
jj.omegastar_x <- jj.discard.psi2$omegastar_x
jj.sigma2squared <- jj.discard.psi2$sigma2squared
jj.omega_t.upper <- chol(jj.omega_t)
jj.omega_t.lower <- t(jj.omega_t.upper)
jj.omega_x.upper <- chol(jj.omega_x)
jj.omega_x.lower <- t(jj.omega_x.upper)
jj.a <- solve(solve(jj.V_theta) + 2 * jj.omega_t, solve(jj.V_theta,
jj.mean))
jj.b <- t(2 * solve(solve(jj.V_theta) + 2 * jj.omega_t) %*%
jj.omega_t)
jj.c <- jj.sigma1squared/sqrt(det(diag(nrow = nrow(jj.V_theta)) +
2 * jj.V_theta %*% jj.omega_t))
names(jj.c) <- "ht.fun.precalc"
jj.A <- solve(jj.V_theta + solve(jj.omega_t)/4)
jj.A.upper <- chol(jj.A)
jj.A.lower <- t(jj.A.upper)
list(rho = rho, lambda = lambda, psi1 = psi1, psi1.apriori = psi1.apriori,
psi2 = psi2, psi2.apriori = psi2.apriori, theta.apriori = theta.apriori,
omega_x = jj.omega_x, omega_t = jj.omega_t,
omegastar_x = jj.omegastar_x, sigma1squared = jj.sigma1squared,
sigma2squared = jj.sigma2squared, omega_x.upper = jj.omega_x.upper,
omega_x.lower = jj.omega_x.lower, omega_t.upper = jj.omega_t.upper,
omega_t.lower = jj.omega_t.lower, a = jj.a, b = jj.b,
c = jj.c, A = jj.A, A.upper = jj.A.upper, A.lower = jj.A.lower)
}
###################################################
### code chunk number 12: make_jj.psi1
###################################################
jj.psi1 <- 1:5
names(jj.psi1) <- c("x", "y", "A","B", "s1sq")
jj.psi2 <- 1:3
names(jj.psi2) <- c("x", "y", "s1sq")
jj.mean1 <- rep(1,5)
names(jj.mean1) <- names(jj.psi1)
jj.sigma1 <- diag(c(1.1, 1.1, 1.2, 1.3, 1.1))
rownames(jj.sigma1) <- names(jj.psi1)
colnames(jj.sigma1) <- names(jj.psi1)
jj.mean2 <- c(1,0.1,rep(1.1,3))
names(jj.mean2) <- c("rho","lambda",names(jj.psi2))
jj.sigma2 <- diag(c(1,0.2,1.1, 1.1, 1.2))/10
rownames(jj.sigma2) <- names(jj.mean2)
colnames(jj.sigma2) <- names(jj.mean2)
jj.mean.th <- 1:2
names(jj.mean.th) <- c("A","B")
jj.sigma.th <- diag(c(1.5,1.7))
rownames(jj.sigma.th) <- names(jj.mean.th)
colnames(jj.sigma.th) <- names(jj.mean.th)
###################################################
### code chunk number 13: make_phi.int
###################################################
phi.int <-
phi.fun.int(rho=1,
lambda=1,
psi1 = jj.psi1,
psi2 = jj.psi2,
psi1.apriori = list(mean=jj.mean1,sigma=jj.sigma1),
psi2.apriori = list(mean=jj.mean2,sigma=jj.sigma2),
theta.apriori = list(mean=jj.mean.th,sigma=jj.sigma.th)
)
###################################################
### code chunk number 14: call_phi.change
###################################################
phi.int2 <- phi.change(old.phi=phi.int, phi.fun=phi.fun.int, rho=3)
print(phi.int2$rho)
###################################################
### code chunk number 15: examine_E.theta.toy
###################################################
E.theta.toy
###################################################
### code chunk number 16: examine_h1.int
###################################################
h1.int
###################################################
### code chunk number 17: make_E.theta.int
###################################################
E.theta.int <- function(D2 = NULL, H1 = NULL, x1 = NULL, x2 = NULL, phi, give.mean = TRUE)
{
if (give.mean) {
m_theta <- phi$theta.apriori$mean
return(H1(D1.fun(D2, t.vec = m_theta)))
}
else {
out <- matrix(0, 5, 5)
out[3,3] <- phi$theta.apriori$sigma[1,1]
out[3,4] <- phi$theta.apriori$sigma[1,1]*x1[1]
out[4,3] <- phi$theta.apriori$sigma[1,1]*x2[1]
out[4,4] <- phi$theta.apriori$sigma[1,1]*x1[1]*x2[1]
out[5,5] <- phi$theta.apriori$sigma[2,2]*x1[2]*x2[2]
return(out)
}
}
###################################################
### code chunk number 18: make_Edash.theta.toy
###################################################
Edash.theta.int <- Edash.theta.toy
###################################################
### code chunk number 19: make_theta.TRUE
###################################################
theta.TRUE <- c(1,1)
###################################################
### code chunk number 20: make_beta1.TRUE
###################################################
beta1.TRUE <- c(0,0,0,1,1)
psi1.TRUE <- c(4,4,4,4,0.5)
###################################################
### code chunk number 21: make_two.designs
###################################################
two.designs <- rbind(D1.int,D1.fun(x.star=D2.int,t.vec=theta.TRUE))
###################################################
### code chunk number 22: examine_two.designs
###################################################
two.designs[c(1:3,(n1+n2):(n1+n2-2)),]
###################################################
### code chunk number 23: make_y.int
###################################################
jj.mean <- H1.int(two.designs) %*% beta1.TRUE
jj.sigma <- psi1.TRUE[5]*corr.matrix(two.designs, scales = psi1.TRUE[1:4])
code.and.obs <- as.vector(rmvnorm(n=1,mean=jj.mean,sigma=jj.sigma))
y.int <- code.and.obs[1:n1]
z.int <- code.and.obs[(n1+1):(n1+n2)]
names(y.int) <- rownames(D1.int)
head(y.int)
###################################################
### code chunk number 24: make_beta2.TRIE
###################################################
beta2.TRUE <- c(1,1)
psi2.TRUE <- c(3,3,0.6)
###################################################
### code chunk number 25: make_model.inadequacy
###################################################
jj.mean <- drop(H2.int(D2.int) %*% beta2.TRUE)
jj.sigma <- corr.matrix(D2.int, scales=psi2.TRUE[1:2])*psi2.TRUE[3]
model.inadequacy <- rmvnorm(n=1, mean=jj.mean,sigma=jj.sigma)
z.int <- as.vector(z.int + model.inadequacy)
names(z.int) <- rownames(D2.int)
###################################################
### code chunk number 26: add_error
###################################################
lambda.TRUE <- 0.00
jj.obs.error <- rnorm(n2)*lambda.TRUE
z.int <- z.int + jj.obs.error
head(z.int)
###################################################
### code chunk number 27: make_d.int
###################################################
d.int <- c(y.int , z.int)
###################################################
### code chunk number 28: make_psi1
###################################################
phi.true <- phi.change(phi.fun=phi.fun.int, old.phi=phi.int,
psi1=psi1.TRUE, psi2=psi2.TRUE,lambda=0.1,rho=1)
###################################################
### code chunk number 29: call_betahat.fun.koh
###################################################
betahat.fun.koh(theta=theta.TRUE, d=d.int, D1=D1.int, D2=D2.int, H1=H1.int, H2=H2.int, phi=phi.true)
###################################################
### code chunk number 30: call_betahat.fun.koh_otherparams
###################################################
betahat.fun.koh(theta=c(-5,5), d=d.int, D1=D1.int, D2=D2.int, H1=H1.int, H2=H2.int, phi=phi.true)
###################################################
### code chunk number 31: call_p.eqn8.supp
###################################################
p.eqn8.supp(theta=c(1,1), D1=D1.int, D2=D2.int, H1=H1.int, H2=H2.int, d=d.int, phi=phi.true)
###################################################
### code chunk number 32: call_p.eqn8.supp_difftheta
###################################################
p.eqn8.supp(theta=c(5,-6), D1=D1.int, D2=D2.int, H1=H1.int, H2=H2.int, d=d.int, phi=phi.true)
###################################################
### code chunk number 33: call_phi.stage1
###################################################
phi.stage1 <- stage1(D1=D1.int, y=y.int, H1=H1.int, maxit=10,
method="SANN", trace=0, do.print=FALSE, phi.fun=phi.fun.int,
phi=phi.int)
###################################################
### code chunk number 34: examine_stage1
###################################################
phi.stage1$psi1
###################################################
### code chunk number 35: examine_stage1_cont
###################################################
phi.stage1$sigma1squared
phi.stage1$omega_x
phi.stage1$omega_t
###################################################
### code chunk number 36: call_stage2
###################################################
use1 <- 1:10
use2 <- 1:11
phi.stage2 <- stage2(D1=D1.int[use1,], D2=D2.int[use2,], H1=H1.int, H2=H2.int,
y=y.int[use1], z=z.int[use2], extractor=extractor.int,
phi.fun=phi.fun.int, E.theta=E.theta.int, Edash.theta=Edash.theta.int,
maxit=1, method="SANN", phi=phi.stage1)
###################################################
### code chunk number 37: examine_stage2
###################################################
phi.stage2$rho
phi.stage2$lambda
phi.stage2$sigma2squared
phi.stage2$psi2
###################################################
### code chunk number 38: make_Xdist
###################################################
jj.xdist.mean <- rep(0.5,2)
names(jj.xdist.mean) <- c("x","y")
jj.xdist.var <- 0.05+diag(c(0.1,0.1))
rownames(jj.xdist.var) <- c("x","y")
colnames(jj.xdist.var) <- c("x","y")
X.dist.int <- list(mean=jj.xdist.mean,var=jj.xdist.var)
###################################################
### code chunk number 39: make_hbar.fun.int
###################################################
hbar.fun.int <-
function (theta, X.dist, phi)
{
if (is.vector(theta)) {
theta <- t(theta)
}
first.bit <- phi$rho * H1.int(D1.fun(X.dist$mean, theta))
first.bit[,5] <- first.bit[,5] + theta[,2]*X.dist$var[2,2]
second.bit <- H2.int(X.dist$mean)
jj.names <- colnames(second.bit)
second.bit <- kronecker(second.bit, rep(1, nrow(first.bit)))
colnames(second.bit) <- jj.names
return(t(cbind(first.bit, second.bit)))
}
###################################################
### code chunk number 40: call_EK.eqn10.theta
###################################################
if(do_from_scratch){
jj <- EK.eqn10.supp(X.dist=X.dist.int, D1=D1.int, D2=D2.int,
H1=H1.int, H2=H2.int, d=d.int, hbar.fun=hbar.fun.int,
lower.theta=c(-3,-3), upper.theta=c(3,3),
extractor=extractor.int,
phi=phi.stage2,eps=0.8)
} else {
jj <- 1.95633284138600
}
###################################################
### code chunk number 41: printjj
###################################################
jj
###################################################
### code chunk number 42: defineMHwrapper
###################################################
p <- function(theta){
p.eqn8.supp(theta=theta, D1=D1.int, D2=D2.int, H1=H1.int, H2=H2.int, d=d.int, phi=phi.true)
}
p(c(1,1))
p(c(10,10))
###################################################
### code chunk number 43: sample.from.MH
###################################################
ss <- MH(n=10,start=c(1,1),sigma=diag(2),pi=p)
###################################################
### code chunk number 44: scatterplot
###################################################
plot(jitter(ss),xlab="A",ylab="B",main="Sample from posterior: try with more points")
###################################################
### code chunk number 45: SetTheBib
###################################################
bib <- system.file( "doc", "uncertainty.bib", package = "emulator" )
bib <- sub('.bib$','',bib)
###################################################
### code chunk number 46: usethebib
###################################################
cat( "\\bibliography{",bib,"}\n",sep='')
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.