3

I'm having an issue where my specific inputs generate an error.

library(kriging)

x <- c(0.65,0.45,0.25,0.65,0.45,0.25,0.55,0.4,0.25,0.55,0.4,0.25,0.6,0.45,0.25,0.6,0.45,0.25,0.5,0.4,0.25,0.5,0.4,0.25,0.55,0.4,0.25,0.55,0.4,0.25,0.5,0.35,0.25,0.5,0.35,0.25)

y <- c(0.25,0.45,0.65,0.2,0.4,0.6,0.25,0.4,0.55,0.2,0.35,0.5,0.25,0.4,0.6,0.2,0.35,0.55,0.25,0.35,0.5,0.2,0.3,0.45,0.25,0.4,0.55,0.2,0.35,0.5,0.2,0.35,0.45,0.15,0.3,0.4)

r <- c(241.5,236.8333333,229.875,242,235.5,231.3333333,238,236.875,225.75,238.5,233.25,228.5,24,0.5,237.1666667,229.5,241.3333333,236.8333333,227.75,237.625,233.7,228.3333333,236.8,235,229.8333333,238.1,234.6,228.6666667,237.375,235.1,228.6666667,236.4,231.6666667,227.3,236.5,232.625,227.3571429)


x <- data.matrix(x)
x <- sweep(x,1,100,"*")

y <- data.matrix(y)
y <- sweep(y,1,100,"*")

krig <- kriging(x,y,r,lag=3)

Error in solve.default(matrix(A, n + 1, n + 1)) :Lapack routine dgesv: system is exactly singular: U[26,26] = 0

I honestly don't know what's wrong. My formatting seems to be well aligned with example code I've seen.

Andrew
  • 7,201
  • 5
  • 25
  • 34
  • 1
    Don't know the details of that package, but my guess is you have duplicate data points which is causing the kernel matrix to be singular, and the package isn't equipped to deal this this situation. – mrip Mar 16 '15 at 20:30
  • Your use of `data.matrix` and `sweep` seems a bit obscure. `x = x * 100` yields the same result... – Paul Hiemstra Mar 16 '15 at 20:54

1 Answers1

8

The problem is, as @mrip mentioned, that you have duplicate observations. In this context, points with exactly the same x and y coordinates. The following code reveals the duplicates:

coor = cbind(x, y)
duplicated(coor)
# [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#[13] FALSE  TRUE  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE
#[25]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE
coor[duplicated(coor),]
#     [,1] [,2]
#[1,]   45   40
#[2,]   25   60
#[3,]   25   55
#[4,]   40   35
#[5,]   25   50
#[6,]   55   25
#[7,]   40   40
#[8,]   25   55
#[9,]   55   20
#[10,]   40   35
#[11,]   25   50
#[12,]   50   20
#[13,]   25   45

This leads to an uninvertible covariance matrix in the kriging equations, which in turn leads to the Lapack error from dgesv. The solution is to remove the duplicates from the dataset.

The following dataset works:

x = runif(100)
y = runif(100)
z = runif(100)
krig = kriging(x, y, z, lag = 3)
image(krig)

enter image description here

Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149