# # Calibration Example # For STAT8810 # Fall, 2017 # M.T. Pratola ########################################################################################################################## # The data, c/o Santner & Leatherman # # The true physical data is generated as # # yf = 0.1x^2-x+0.4+e # # where e ~ N(0,0.04) # # The simulator outputs are generated as # # yc = theta1*x + theta2 # # In both, x is in [-5,5]. # # Note that this setup means the simulator is biased since it does not include # the quadratic component of the true real-world process. # So we expect the discrepancy to be non-zero. ########################################################################################################################## design=matrix(c( -5.0000,-1.4286,-0.7143, -4.2857, 3.5714,-2.1429, -3.5714,-2.1429, 4.2857, -2.8571, 2.8571, 2.8571, -2.1429,-4.2857,-3.5714, -1.4286, 0.7143,-5.0000, -0.7143,-5.0000, 1.4286, 0.0000, 0.0000, 0.0000, 0.7143, 5.0000,-1.4286, 1.4286,-0.7143, 5.0000, 2.1429, 4.2857, 3.5714, 2.8571,-2.8571,-2.8571, 3.5714, 2.1429,-4.2857, 4.2857,-3.5714, 2.1429, 5.0000, 1.4286, 0.7143),ncol=3,byrow=TRUE) yc=c(6.4286,-17.4490,11.9388,-5.3061,5.6122,-6.0204,5.0000,0,2.1429,3.9796,12.7551,-11.0204,3.3673,-13.1633,7.8571) th.init=0 #always d.th=matrix(c( -5.0000,th.init,th.init, -2.5000,th.init,th.init, 0.0000,th.init,th.init, 2.5000,th.init,th.init, 5.0000,th.init,th.init),ncol=3,byrow=TRUE) yf=c(7.6300,4.1320,0.5451,-1.4876,-1.9571) my=mean(yc) sdy=sd(yc) yc=(yc-my)/sdy yf=(yf-my)/sdy r1=range(design[,1]) r2=range(design[,2]) r3=range(design[,3]) design[,1]=(design[,1]-r1[1])/(r1[2]-r1[1]) design[,2]=(design[,2]-r2[1])/(r2[2]-r2[1]) design[,3]=(design[,3]-r3[1])/(r3[2]-r3[1]) d.th[,1]=(d.th[,1]-r1[1])/(r1[2]-r1[1]) nf=length(yf) nc=length(yc) y=c(yf,yc) design.all=rbind(d.th,design) l1=list(m1=outer(design.all[,1],design.all[,1],"-")) l2=list(m2=outer(design.all[,2],design.all[,2],"-")) l3=list(m3=outer(design.all[,3],design.all[,3],"-")) l.dez=list(l1=l1,l2=l2,l3=l3) ncdims=2 nxdims=1 ninfo=list(npred=0,nf=nf,nc=nc,ncdims=ncdims) # finally set up the prediction grid: npred=30 design.pred=as.matrix(expand.grid(seq(0,1,length=npred),th.init,th.init)) dpred.all=rbind(design.pred,d.th,design) l1=list(m1=outer(dpred.all[,1],dpred.all[,1],"-")) l2=list(m2=outer(dpred.all[,2],dpred.all[,2],"-")) l3=list(m3=outer(dpred.all[,3],dpred.all[,3],"-")) l.pred=list(l1=l1,l2=l2,l3=l3) ninfo.pred=list(npred=npred,nf=nf,nc=nc,ncdims=ncdims) ########################################################################################################################## # Now fit the calibration model ########################################################################################################################## source("dace.sim.r") source("calibration.r") pi=list(af=80,bf=4,az=5,bz=5,adel=12,bdel=2,rhoa=rep(5,nxdims),rhob=rep(1,nxdims),psia=rep(5,ncdims),psib=rep(1,ncdims),rdla=rep(5,nxdims),rdlb=rep(1,nxdims)) mh=list(rr=.03,rp=.1,re=.1,rf=50,rz=.5,rl=.7,rth=.15) set.seed(88) # just to replicate my run fit=calibrate(y,l.dez,ninfo,20000,pi,mh,last=10000) fitp=predict(fit,l.pred,ninfo.pred) ########################################################################################################################## # Finally make some plots. ########################################################################################################################## bias=fitp$z.pred-fitp$eta.pred bias.mu=apply(bias,2,mean) bias.sd=apply(bias,2,sd) #plot priors pdf("prior.caleg.pdf",width=10,height=10) par(mfrow=c(2,2)) plot(density(rgamma(1e5,shape=pi$af,rate=pi$bf)),main="",xlab="lambdaf prior") plot(density(rgamma(1e5,shape=pi$az,rate=pi$bz)),main="",xlab="lambdaz prior") plot(density(rgamma(1e5,shape=pi$adel,rate=pi$bdel)),main="",xlab="lambdadel prior") plot(density(rbeta(1e5,pi$rhoa[1],pi$rhob[1])),main="",xlab="corr. param prior") dev.off() pdf("post.caleg.pdf",width=10,height=10) par(mfrow=c(3,3)) plot(fit$lambdaf,type='l',xlab="lambdaf") plot(fit$lambdaz,type='l',xlab="lambdaz") plot(fit$lambdadel,type='l',xlab="lambdadel") plot(fit$rhoz,type='l',xlab="rhoz",ylim=c(0,1)) plot(fit$psiz[,1],type='l',xlab="psiz.1",ylim=c(0,1)) plot(fit$psiz[,2],type='l',xlab="psiz.2",ylim=c(0,1)) plot(fit$rhodel,type='l',xlab="rhodel",ylim=c(0,1)) plot(fit$theta[,1]*(r2[2]-r2[1])+r2[1],type='l',xlab="theta.1",ylim=c(-5,5),ylab="") plot(fit$theta[,2]*(r3[2]-r3[1])+r3[1],type='l',xlab="theta.2",ylim=c(-5,5),ylab="") dev.off() pdf("post.theta.caleg.pdf",width=10,height=10) par(mfrow=c(1,2)) plot(density(fit$theta[,1]*(r2[2]-r2[1])+r2[1]),xlab="theta.1",xlim=c(-5,5),main="") abline(v=-1,lty=2) plot(density(fit$theta[,2]*(r3[2]-r3[1])+r3[1]),xlab="theta.2",xlim=c(-5,5),main="") abline(v=0.4,lty=2) dev.off() pdf("post.pred.caleg.pdf") par(mfrow=c(1,1)) plot(d.th[,1],yf,pch=20,col="red",ylim=c(-2,2),xlab="x",ylab="response") points(design[,1],yc,pch=20,col="grey") lines(design.pred[,1],fitp$z.mu,col="green") lines(design.pred[,1],fitp$z.mu-1.96*fitp$z.sd,col="green",lty=2) lines(design.pred[,1],fitp$z.mu+1.96*fitp$z.sd,col="green",lty=2) lines(design.pred[,1],fitp$eta.mu,col="blue") lines(design.pred[,1],fitp$eta.mu-1.96*fitp$eta.sd,col="blue",lty=2) lines(design.pred[,1],fitp$eta.mu+1.96*fitp$eta.sd,col="blue",lty=2) lines(design.pred[,1],bias.mu,col="orange") lines(design.pred[,1],bias.mu-1.96*bias.sd,col="orange",lty=2) lines(design.pred[,1],bias.mu+1.96*bias.sd,col="orange",lty=2) lines(design.pred[,1],0.1*design.pred[,1]^2,lty=4) #true discrepancy dev.off() cortb=cor(fit$theta[,2],apply(bias,1,mean)) # -0.88 highly correlated pdf("post.ident.caleg.pdf",width=10,height=5) par(mfrow=c(1,1)) plot(fit$theta[,2],apply(bias,1,mean),xlab="theta.2",ylab="E[bias]") #plot the theta.2 draws vs the mean (over xpred) bias to show the correlation dev.off()