I'm plotting Matern correlation function with various values of the shape parameter (0.5, 1.5, and 2.5) and scale parameter (0.25, 0.16, and 0.13, respectively). The values are motivated by an example in Diggle & Ribeiro (2007) Model-based Geostatistics book, which is related to the geoR package. I tried to plot the correlation functions based on their geoR example, and then tried to make the same plots using RandomFields package. The plots suppose to be same, but their are not (see below). I checked the parametrization, and it seems correct.
library(geoR)
x <- seq(0, 1.5, l=200)
m1.geoR <- cov.spatial(x, cov.model = "mat", kappa=0.5, cov.pars = c(1, 0.25))
m2.geoR <- cov.spatial(x, cov.model = "mat", kappa=1.5, cov.pars = c(1, 0.16))
m3.geoR <- cov.spatial(x, cov.model = "mat", kappa=2.5, cov.pars = c(1, 0.13))
plot(x, m1.geoR, "l")
lines(x, m2.geoR, lty=2)
lines(x, m3.geoR, lty=3)
library(RandomFields)
m1.RF <- RMmatern(nu=0.5, scale=0.25, var = 1)
m2.RF <- RMmatern(nu=1.5, scale=0.16, var = 1)
m3.RF <- RMmatern(nu=2.5, scale=0.13, var = 1)
plot(m1.RF)
lines(m2.RF, lty=2)
lines(m3.RF, lty=3)
Plot generated using RandomFields
I also tried to use my own function for the Matern model, and the resulting plot corresponds to the plot based on geoR. So it seems there is something wrong in the RandomFields implementation. Or am I missing something? I was not able to find the actual implementation used in the RandomFields package. Any explanation would be greatly appreciated.
matern <- function(x, nu, scale){
ifelse(x==0, 1, (2^(1-nu))*(gamma(nu)^(-1))*(x/scale)^nu*besselK((x/scale), nu))
}
m1.matern <- matern(x, nu=0.5, scale=0.25)
m2.matern <- matern(x, nu=1.5, scale=0.16)
m3.matern <- matern(x, nu=2.5, scale=0.13)
plot(x, m1.matern, "l")
lines(x, m2.matern, lty=2)
lines(x, m3.matern, lty=3)
