require(geoR)
set.seed(285)
##
## Further examples for the usage of krige.bayes()
## -----------------------------------------------
##
## Before reading this see documentation for krige.bayes
## > help(krige.bayes)
##
## WARNING: THESE CALLS CAN BE TIME AND SPACING DEMMANDING

## Simulating data
ex.data <- grf(50, cov.pars=c(10, .25))

##
## Basic usage
##

## a basic and simple call to the function
ex.post <- krige.bayes(ex.data)
ex.post
names(ex.post)

## different input and output options
ex1 <- krige.bayes(ex.data, prior = list(phi.prior = "fixed", phi = 0.3))
ex1 <- krige.bayes(ex.data, model = list(cov.model="spherical"))
ex1 <- krige.bayes(ex.data, output = list(n.posterior = 100))

## now performing prediction
ex.grid <- as.matrix(expand.grid(seq(0,1,l=6), seq(0,1,l=6)))
ex.bayes <- krige.bayes(ex.data, loc=ex.grid, prior =
                        prior.control(phi.discrete=seq(0, 2, l=3),
                                      tausq.rel.discrete=seq(0, 2, l=3)),
                        output=output.control(n.post=100))
names(ex.bayes)

## some graphical output ...
## ... for the posterior ...
plot(ex.data)
lines(ex.bayes, sum = mean) 
plot(ex.bayes)
lines(ex.bayes, summ="median", lty=2, post="par")
lines(ex.bayes, summ="mean", lwd=2, lty=2, post="par")

## ... and for the predictive
op <- par(no.readonly = TRUE)
par(mfrow=c(2,2))
par(mar=c(3,3,1,1))
par(mgp = c(2,1,0))
image(ex.bayes, main="predicted values")
image(ex.bayes, val="variance", main="prediction variance")
image(ex.bayes, val= "simulation", number.col=1,
      main="a simulation from the \npredictive distribution")
image(ex.bayes, val= "simulation", number.col=2,
      main="another simulation from \nthe predictive distribution")

## Now including tausq (the nugget variance term)
PC <- prior.control(phi.prior = c(.1, .2, .3, .2, .1, .1), phi.disc=seq(0.1, 0.6, l=6),
                    tausq.rel.prior=c(.1, .4, .3, .2), tausq.rel.discrete=c(0,.1,.2,.3))
ex.user <- krige.bayes(ex.data, prior = PC)  
ex.user

##
## Example for "sequential analysis"
##
## Simulating data at 2 different "times"
ap1 <- grf(50, cov.pars=c(1, .3))
ap2 <- grf(70, cov.pars=c(1, .3))
## A initial "usual" analysis
ap1.kb <- krige.bayes(ap1)
ap1.kb
## and now using the previously obtained posterior as prior for next call
ap2.kb <- krige.bayes(ap2, prior=post2prior(ap1.kb))
ap2.kb

##
## Another example with a "user defined" prior
##
PC <- prior.control(phi.prior=c(.2,.3,.2,.1,.1,.1), phi.discrete = seq(0,.5,l=6))
##ap3.kb <- krige.bayes(ap1, prior = PC)
ap3.kb <- krige.bayes(ap1, prior = PC)
ap3.kb
##
ap4.kb <- krige.bayes(ap2, prior=post2prior(ap3.kb))
ap4.kb

##
## Now including tausq
##
PC <- prior.control(tausq.rel.prior = "uni", tausq.rel.discrete = seq(0, .5, l=6))
ap5.kb <- krige.bayes(ap1)
ap5.kb
## using the previous posterior as prior for next call
ap6.kb <- krige.bayes(ap2, prior=post2prior(ap5.kb))
ap6.kb
##
##
##
PC <- prior.control(phi.prior=c(.2,.3,.2,.1,.1,.1), phi.discrete = seq(0,.5,l=6),
                    tausq.rel.prior=c(.3,.4,.3), tausq.rel.discrete = c(0, .1, .2)) 
ap7.kb <- krige.bayes(ap1, prior = PC)
ap7.kb

ap8.kb <- krige.bayes(ap2, prior=post2prior(ap7.kb))
ap8.kb

##
## An example with a trend term in the model and different priors
##
data(s100)
prior2.b9 <- prior.control(beta.prior = "normal", beta = c(0,0,0),
                           beta.var = cbind(c(2,1.5,0),c(1.5,1.8,.2),c(0,0.2,1.5)),
                           phi.prior = "exponential", phi = 2.5, phi.discrete = c(2.5,3),
                           sigmasq.prior = "sc.inv.chisq", df.sigmasq = 5, sigmasq = 0.5)
ap <- krige.bayes(s100, prior=prior2.b9, model=model.control(trend.d = "1st"))
ap

prior2.b9 <- prior.control(beta.prior = "normal", beta = c(0,0,0),
                           beta.var = cbind(c(2,1.5,0),c(1.5,1.8,0.5),c(0,0.5,1.5)),
                           phi.prior = "fixed", phi = 2.5, sigmasq.prior = "sc.inv.chisq",
                           df.sigmasq = 5, sigmasq = 0.5)
ap <- krige.bayes(s100, prior=prior2.b9, model=model.control(trend.d="1st"))
ap

prior2.b9 <- prior.control(beta.prior = "normal", beta = 0,beta.var =1,
                           phi.prior = "fixed", phi = 2.5,sigmasq.prior ="sc.inv.chisq",
                           df.sigmasq = 5, sigmasq = 0.5)
ap <- krige.bayes(s100, prior=prior2.b9)
ap

par(op)