rmatlabstatisticsspatstatstochastic-process

Comparing estimates of Ripley's K function using Matlab and R


I am using the following Matlab code to estimate Ripley's K function.

a = 0;
b = 50;
C_x = a + (b-a).*rand(100,1);
C_y = a + (b-a).*rand(100,1);

locs = zeros(length(C_x),2);
locs(:,1) = C_x;
locs(:,2) = C_y;

dist = a:1:b;
K_t = RipleysK(locs, dist, [a, b, a, b]);
plot(dist,K_t);
title('$\hat{K}$ function','Interpreter','latex','FontSize',14);
xlabel('$r$','Interpreter','latex','FontSize',14);
ylabel('$\hat{K}(r)$','Interpreter','latex','FontSize',14);
csvwrite('test.csv',locs);

"RipleysK" function can be found at: http://www.colorado.edu/geography/class_homepages/geog_4023_s07/labs/lab10/RipleysK.m

In comparison, I am using the following R script.

mydata <- read.csv("test.csv",header=F)
mydata
w <- owin(xrange = c(0,50),yrange = c(0,50))
pp <- as.ppp(mydata,w)
K <- Kest(pp,correction = "none")
plot(K)

Matlab estimates K for specified values of r (i.e., dist) whereas R script does not (estimates till r = 12.5).

Can somebody comment please? Which one is correct? Can we specify r values in R script?

Thanks


Solution

  • I have not looked the matlab code, but I'm very certain that Kest in the spatstat package calculates the right thing. From your own comment it also appears that there is no conflict in this -- they just do two different things: The matlab code evaluates an estimate of the K-function for one specific value of r while Kest does it for a range of r-values (and it also calculates a collection of different estimators). By default it calculates the estimates for 513 values of r, and the easiest way to see them all is to covert to a data.frame (continuing your code):

    Kdf <- as.data.frame(K)
    head(Kdf)
    

    Alternatively you can convert the function value (fv) object K to a normal R function and evaluate that at the relevant value of r (this automatically chooses the "best" estimator among the available ones (typically Ripley's isotrotically corrected estimator)):

    Kf <- as.function(K)
    Kf(12.5)