if (require(INLA)) {
#Generate mesh and SPDE model
n.lattice = 20 #increase for more interesting, but slower, examples
x=seq(from=0,to=10,length.out=n.lattice)
lattice=inla.mesh.lattice(x=x,y=x)
mesh=inla.mesh.create(lattice=lattice, extend=FALSE, refine=FALSE)
spde <- inla.spde2.matern(mesh, alpha=2)
#Generate an artificial sample
sigma2.e = 0.01
n.obs=1000
obs.loc = cbind(runif(n.obs)*diff(range(x))+min(x),
runif(n.obs)*diff(range(x))+min(x))
Q = inla.spde2.precision(spde, theta=c(log(sqrt(0.5)), log(sqrt(1))))
x = inla.qsample(Q=Q)
A = inla.spde.make.A(mesh=mesh,loc=obs.loc)
Y = as.vector(A %*% x + rnorm(n.obs)*sqrt(sigma2.e))
## Calculate posterior
Q.post = (Q + (t(A)%*%A)/sigma2.e)
mu.post = as.vector(solve(Q.post,(t(A)%*%Y)/sigma2.e))
## Calculate continuous contours
tric = tricontour(mesh, z=mu.post, levels=as.vector(quantile(x,c(0.25,0.75))))
}Run the code above in your browser using DataLab