rMatern <- function(n, coords, kappa, variance, nu=1) { m <- as.matrix(dist(coords)) m <- exp((1-nu)*log(2) + nu*log(kappa*m)- lgamma(nu))*besselK(m*kappa, nu) diag(m) <- 1 return(drop(crossprod(chol(variance*m), matrix(rnorm(nrow(coords)*n), ncol=n)))) } rspde <- function(coords, kappa, variance=1, alpha=2, n=1, mesh, verbose=FALSE, seed, return.attributes=FALSE) { t0 <- Sys.time() theta <- c(-0.5*log(4*pi*variance*kappa^2), log(kappa)) if (verbose) cat('theta =', theta, '\n') if (missing(mesh)) { mesh.pars <- c(0.5, 1, 0.1, 0.5, 1)*sqrt(alpha-ncol(coords)/2)/kappa if (verbose) cat('mesh.pars =', mesh.pars, '\n') attributes <- list( mesh=inla.mesh.2d(, coords[chull(coords), ], max.edge=mesh.pars[1:2], cutoff=mesh.pars[3], offset=mesh.pars[4:5])) if (verbose) cat('n.mesh =', attributes$mesh$n, '\n') } else attributes <- list(mesh=mesh) attributes$spde <- inla.spde2.matern(attributes$mesh, alpha=alpha) attributes$Q <- inla.spde2.precision(attributes$spde, theta=theta) attributes$A <- inla.mesh.project(mesh=attributes$mesh, loc=coords)$A if (n==1) result <- drop(attributes$A%*%inla.qsample( Q=attributes$Q, constr=attributes$spde$f$extraconstr)) t1 <- Sys.time() result <- inla.qsample(n, attributes$Q, seed=ifelse(missing(seed), 0, seed), constr=attributes$spde$f$extraconstr) if (nrow(result)