2

I was wondering if any one could explain to me how the geoR package calculates the covariance function? I mean how you would do it by hand?

library(geoR)
#suppose I have the following coordinates
X = c(60,30,20,40)
Y = c(50,20,50,50)
my_coordinates = cbind(X,Y)
print(my_coordinates)

#computing covariance
my_cov= varcov.spatial(my_coordinates,cov.model="exp", cov.pars=c(0.2,25))
print(my_cov)

And you get:

         [,1]       [,2]       [,3]       [,4]
[1,] 0.20000000 0.03664442 0.04037930 0.08986579
[2,] 0.03664442 0.20000000 0.05645288 0.05645288
[3,] 0.04037930 0.05645288 0.20000000 0.08986579
[4,] 0.08986579 0.05645288 0.08986579 0.20000000

However, one might want to do it in Matlab as well.

www
  • 38,575
  • 12
  • 48
  • 84
ToNoY
  • 1,358
  • 2
  • 22
  • 43
  • Have you looked at the source? You can do so just by typing the function name at the R command prompt. The function includes several variations according to the supplied args, but the individual code paths don't look particularly obscure. – walkytalky Mar 19 '13 at 22:27
  • As stated below, your distance matrix `h` does not correspond to your data. – mnel Mar 20 '13 at 02:23

1 Answers1

3

The best way to find out how a package or function does something is to look at the source code. This is one of the awesome things about open source projects, you can do this.

try typing varcov.spatial or searching through the unpacked package tar ball for the function definition

To calculate the covariance (which is dependant on the distance between points), you need to calculate

  • the distance between your points (you really only need the lower triangle, as it will be symmetric
  • The value of the covariance function at each distance
  • form the full symmetric variance covariance matrix from these calculated covariances.

The covariance functions are defined in ?cov.spatial. You can call cov.spatial to calculate these in R (exactly what geoR::varcov.spatial does)

mnel
  • 113,303
  • 27
  • 265
  • 254
  • @ToNoY -- your example distances aren't from your example data. The lower triange of your distance matrix is `as.vector(dist(my_coordinates))` which gives `[1] 42.42641 40.00000 20.00000 31.62278 31.62278 20.00000` – mnel Mar 20 '13 at 01:56
  • @ToNoY -- the default is [euclidian distance](http://en.wikipedia.org/wiki/Euclidean_distance) -- see `?dist`. I don't know where you got your numbers from. – mnel Mar 20 '13 at 02:43