4

I'm trying to fit a variogram to my data, however the spatial points are confined by an irregular polygon. So I'd like to supply a variogram model function with the distance matrix of the points.

I've looked through gstat, geoR and vardiag but I can't figure out how to specify the distance matrix, as each of these requires the coordinates of the points.

So I'm really just looking for a way to supply these distances to a variogram model fitting function, any help is appreciated!

smccain
  • 659
  • 5
  • 19
  • 1
    Looking into the `variog` function in `geoR`, I think it might be as simple as adding another argument which changes the distance matrix. In the function: `u – smccain Mar 24 '15 at 17:37
  • It might be a good idea to modify the function accordingly and e-mail the package maintainer with the updated version. At the very least, you can (and should) post a link to your modified function as an answer. On the other hand, it might be unsatisfying to award the bounty to yourself. – shadowtalker Mar 27 '15 at 18:26
  • I've modified the function incorporating another argument which would specify a distance matrix. However I'm running into an error message `could not find function ".define.bins"`. I can't seem to find this function on the help pages - any idea where this comes from? – smccain Mar 27 '15 at 20:14
  • Try `geoR:::.define.bins`, or look in the source code for `geoR` to see if that function gets defined locally. – shadowtalker Mar 27 '15 at 20:16

1 Answers1

0

I don't know if there is a "proper" way to do this without changing the code, as you suggest in your comment; I haven't looked. But here's a useful hack: just make up points whose distances are the same as the distance matrix, and pass those.

  1. Let your squared distance matrix be $D$, of shape $n \times n$. We're looking for an embedding $\{ x_i \}_{i=1}^n$, with $x \in \mathbb R^d$ for some $d$, such that $D_{ij} = \lVert x_i - x_j \rVert^2$.
  2. Find a corresponding Gram matrix $G_{ij} = x_i^T x_j$. Note that $D_{ij} = \lVert x_i \rVert^2 + \lVert x_j \rVert^2 - 2 G_{ij}$. Defining $H$ to be the centering matrix $I - \mathbf 1 \mathbf 1^T / n$, $G = - \frac12 H D H$ works.
  3. Take the eigendecomposition $G = Q \Lambda Q^T$. All the eigenvalues will be nonnegative, since Gram matrices are positive semidefinite; if an exact embedding into dimension $d$ is possible, there will be at most $d$ nonzero entries. Otherwise, truncating it to the highest $d$ entries gives the best (in a certain sense) embedding into dimension $d$. Let $\Lambda_d$ be the diagonal matrix of the $d$ highest eigenvalue, and $Q_d$ be the corresponding eigenvectors (columns of $Q$); thus (assuming an exact embedding is possible) $G = Q_d \Lambda_d Q_d^T$. But then we can just use the embedding $X = Q_d \Lambda_d^{1/2}$, where since $\Lambda_d$ is diagonal the square root is just elementwise: then $X X^T = Q_d \Lambda_d^{1/2} \Lambda_d^{1/2} Q_d^T = G$.
Danica
  • 21,852
  • 1
  • 59
  • 115
  • 1
    This is quite clever--and without changing the code, is likely the best possible approach if it can be made to work. But almost all variography software accepts points of limited dimensions, typically less than $n=2$ or $n=3$; never more (as far as I am aware). Thus you effectively have to reconstruct the original point configuration (up to a Euclidean motion) from the distance matrix. – whuber Mar 27 '15 at 17:53
  • @whuber Good point; I'm not at all familiar with variograms, as was perhaps obvious. But if we use the eigendecomposition instead, as in my edit, we should be able to get down to the source number of dimensions rather than the sillily inflated dimension from using Cholesky that I had before. – Danica Mar 27 '15 at 18:16