In R, I am trying to calculate a Matérn covariance matrix whose inputs are randomly created distance matrices. However, I end up getting covariance matrices that are not positive definite, which makes little sense.
First I generate a random distance matrix:
nrows <- 100
ncols <- 100
d <- matrix(runif(nrows*ncols, 0, 1), ncols, nrows)
# enforce symmetry:
d[lower.tri(d)] <- t(d)[lower.tri(d)]
diag(d) <- 0.000000000001 # instead of zero to avoid numerical issues
# check that we have a correct distance matrix:
range(d)
Then I follow this Wikipedia formula for Matérn covariances: enter image description here
# Calculate Matérn covariance matrix:
sigma <- 1
v <- 3
p <- 5
term1 <- (2**(1-v))/(gamma(v))
term2 <- (sqrt(2*v)*(abs(d)/p))**v
term3 <- besselK(sqrt(2*v)*(abs(d)/p), nu = v)
m <- (sigma**2)*term1*term2*term3
However, when I check the eigenvalues of the Matérn covariance matrix, I see that it's not a positive definite matrix, as it should be:
g <- eigen(m, only.values=TRUE)
print(min(g$values)) # should be greter than zero
print(sum(g$values<0)/length(g$values)) # thus, this should be zero
I fail to understand what am I doing wrong and how to make sure I generate a positive definite Matérn covariance matrix when passing to it a distance matrix.