INLA spatial parameters test

gistfile1.r

library(INLA)
m = 50
points = matrix(runif(m*2), m, 2)
mesh = inla.mesh.2d(loc=points,cutoff=0.05,offset=c(0.1, 0.4),max.edge=c(0.05, 0.5) )
bnd = inla.nonconvex.hull(points, convex=0.12)
mesh = inla.mesh.2d(boundary=bnd,cutoff=0.05,offset = c(1, 0.5),max.edge=c(0.1, 0.5) )
A = inla.spde.make.A(mesh, loc=points)
sigma0 = 1 ## Field std.dev. for theta=0
size = min(c(diff(range(mesh$loc[,1])), diff(range(mesh$loc[,2]))))
range0 = size/5 ## A fifth of the approximate domain width.
kappa0 = sqrt(8)/range0
tau0 = 1/(sqrt(4*pi)*kappa0*sigma0)
spde = inla.spde2.matern(mesh,
                         B.tau=cbind(log(tau0), -1, +1),
                         B.kappa=cbind(log(kappa0), 0, -1),
                         theta.prior.mean=c(0,0),
                         theta.prior.prec=c(0.1, 1) )
mesh.index = inla.spde.make.index(name="field", n.spde=spde$n.spde)

sigma.true = 3 * sigma0
range.true = range0
theta.true = log(c(sigma.true/sigma0, range.true/range0))
Q = inla.spde.precision(spde, theta=theta.true)

x = inla.qsample(n=1, Q)
y = as.vector(A %*% x)
stack = inla.stack(data=list(y=y), A=list(A), effects=list(mesh.index))
result = inla(y ~ -1 + f(field, model=spde), family="gaussian",
              data=inla.stack.data(stack, spde=spde),
              control.predictor=list(A=inla.stack.A(stack)),
              control.family=list(initial=0, fixed=TRUE))
summary(result)

spde2.result = inla.spde2.result(result, "field", spde)
plot(spde2.result$marginals.range.nominal[[1]]) # Nominal range posterior density
abline(v=range.true)
plot(spde2.result$marginals.variance.nominal[[1]]) # Nomial variance posterior density
abline(v=sigma.true^2)

# Convert thetas to variance and range
c(sigma.est = exp(spde2.result$summary.theta$mean[1]) * sigma0, sqrt(exp(spde2.result$summary.log.variance.nominal$mean)))
c(range.est = exp(spde2.result$summary.theta$mean[2]) * range0, exp(spde2.result$summary.log.range.nominal$mean))