0

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

MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
Atta Ul Mustafa
  • 13
  • 2
  • Have you read [`?Kest`](http://www.inside-r.org/packages/cran/spatstat/docs/Kest)? – MichaelChirico Aug 15 '15 at 21:28
  • @MichaelChirico, yes skimmed through R code. I have one observation that Kest estimates K function for r = n/4 e.g., 50/4 = 12.5 which is same as done in Matlab code. They have similar estimates. Now I just wanted to know that whether we can access K values of Kest in R plot (in a similar way we use data tips in Matlab graph) or we can not do that? – Atta Ul Mustafa Aug 16 '15 at 20:42

1 Answers1

0

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)
Ege Rubak
  • 4,347
  • 1
  • 10
  • 18