#1. Plug-in patial predicton for log-transformed lead concentrations

coords <-lead$coords
poly <- bbox(as.matrix(coords))
h<- 10000 # use a smaller value to get a finer grid of prediction points and a nicer-looking map
grid.pred <- expand.grid(seq(poly[1,1],poly[1,2],h),
                         seq(poly[2,1],poly[2,2],h))

# set spatial variance, exponential correlation range and nugget variance at maximum likelihood estimates
sigmasq.and.phi<-MLE$cov.pars
tausq<-MLE$nugget

# set remaining kriging control arguments to perform ordinary kriging at prediction grid-points
kc<-krige.control(cov.model="exponential",cov.pars=sigmasq.and.phi,nugget=tausq,type.krige="ok")

# perform kriging with 100 predictive simulations (would use more than this in a real application)
oc<-output.control(n.predictive=100)
result<-krige.conv(loglead,locations=grid.pred,krige=kc,output=oc)
names(result)
# [1] "predict"              "krige.var"            "beta.est"            
# [4] "distribution"         "simulations"          "mean.simulations"    
# [7] "variance.simulations" "sim.means"            ".Random.seed"        
#[10] "message"              "call" 

predicted.values<-result$mean.simulations
# set predictions outside Galicia to NA (missing)
library(splancs) # required package
predicted.values[!inout(grid.pred,boundary)]<-NA

# draw map
plot(boundary[,1],boundary[,2],type="l",asp=1)
ux<-unique(grid.pred[,1]); nx<-length(ux)
uy<-unique(grid.pred[,2]); ny<-length(uy)
image.plot(ux,uy,matrix(result$mean.simulations,nx,ny),add=TRUE)
lines(boundary[,1],boundary[,2])

#2. Predictive probability of exceedance map

sims<-result$simulations
p<-apply(1*(sims>log(3)),1,mean)
p[!inout(grid.pred,boundary)]<-NA
plot(boundary[,1],boundary[,2],type="l",asp=1)
image.plot(ux,uy,matrix(p,nx,ny),add=TRUE)
lines(boundary[,1],boundary[,2])

#3.Repeat exercises 1 and 2 using Bayesian method of inference

MC<-model.control(kappa=0.5)
PC<-prior.control(beta.prior="flat",sigmasq.prior="sc.inv.chisq",
sigmasq=0.2,df.sigmasq=4,phi.discrete=10000*(1:5),
tausq.rel.prior="uniform",tausq.rel.discrete=0.005*(1:5))
OC<-output.control(n.posterior=100,n.predictive=100,
    simulations.predictive=T,signal=T,moments=F)
result.bayes<-krige.bayes(geodata=loglead,locations=grid.pred,
    model=MC,prior=PC,output=OC)

plot(result.bayes)
posteriors.bayes<-results.bayes$posterior
posterior.sample<-posteriors.bayes$sample
par(mfrow=c(1,3))
hist(posterior.sample[,2],main="sigmasq")
hist(posterior.sample[,3],main="phi")
plot(posterior.sample[,2],posterior.sample[,3],xlab="sigmasq",ylab="phi")

par(mfrow=c(1,1))
predictions.bayes<-result.bayes$predictive
predicted.values.bayes<-predictions.bayes$mean.simulations
predicted.values.bayes[!inout(grid.pred,boundary)]<-NA
plot(boundary[,1],boundary[,2],type="l",asp=1)
ux<-unique(grid.pred[,1]); nx<-length(ux)
uy<-unique(grid.pred[,2]); ny<-length(uy)
image.plot(ux,uy,matrix(predicted.values.bayes,nx,ny),add=TRUE)
lines(boundary[,1],boundary[,2])

# comparison
plot(predicted.values,predicted.values.bayes,pch=19)

sims.bayes<-result.bayes$predictive$simulations
p.bayes<-apply(1*(sims.bayes>log(3)),1,mean)
p.bayes[!inout(grid.pred,boundary)]<-NA
plot(boundary[,1],boundary[,2],type="l",asp=1)
image.plot(ux,uy,matrix(p.bayes,nx,ny),add=TRUE)
lines(boundary[,1],boundary[,2])

# comparison
plot(p,p.bayes,pch=19)