Set weights for ewcdf {spatstat} [R]

Question

I want to compare a reference distribution d_1 with a sample d_2 drawn proportionally to size w_1 using the Kolmogorov–Smirnov distance.

Given that d_2 is weighted, I was considering accounting for this using the Weighted Empirical Cumulative Distribution Function in R (using ewcdf {spatstat}).

The example below shows that I am probably miss-specifying the weights, because when lenght(d_1) == lenght(d_2) the Kolmogorov–Smirnov is not giving a value of 0.

Can someone help me with this? For clarity, see the reproducible example below.

#loop for testing sample sizes 1:length(d_1)
d_stat <- data.frame(1:1000, rep(NA, 1000))
names(d_stat) <- c("sample_size", "ks_distance")

for (i in 1:1000) {

#reference distribution
d_1 <- rpois(1000, 500)
w_1 <- d_1/sum(d_1)
m_1 <- data.frame(d_1, w_1)

#sample from the reference distribution
m_2 <-m_1[(sample(nrow(m_1), size=i, prob=w_1, replace=F)),]
d_2 <- m_2$d_1
w_2 <- m_2$w_1

#ewcdf for the reference distribution and the sample
f_d_1 <- ewcdf(d_1)
f_d_2 <- ewcdf(d_2, 1/w_2, normalise=F, adjust=1/length(d_2))

#kolmogorov-smirnov distance
d_stat[i,2] <- max(abs(f_d_1(d_2) - f_d_2(d_2)))
}

d_stat[1000,2]

Answer 1

I don’t quite understand what you are trying to do here. Why would you expect ewcdf(d_1) and ewcdf(d_2, w_2, normalise=F) to give the same result for i=1000? The first one is the usual ecdf which jumps at the unique values of the input vector with a jump size determined by the number of times the value is repeated (more ties – larger jumps). The second one jumps at the same unique values with a height determined by the sum of the weights you have provided.

What does give identical results is ewcdf(d_2, w_2) and ewcdf(d_1, w_1), but this is not the same as ewcdf(d_1). To understand why the latter two are different, I would suggest a much smaller handmade example with a couple of ties:

library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: nlme
#> Loading required package: rpart
#> 
#> spatstat 1.60-1.006       (nickname: 'See Above') 
#> For an introduction to spatstat, type 'beginner'
x <- c(1,2,3,3,4)
e <- ewcdf(x)

This is the usual ecdf which jumps with value 1/5 at x=1, 1/5 at x=2, 2*1/5 at x=3 and 1/5 at x=4:

plot(e)

Now you define the weights as:

w <- x/sum(x)
w
#> [1] 0.07692308 0.15384615 0.23076923 0.23076923 0.30769231

Thus the ewcdf will jump with value 1/13 at x=1, 2/13 at x=2, 2*3/13 at x=3 and 4/13 at x=4 (with the usual ecdf overlayed in red):

plot(ewcdf(x, w, normalise = FALSE), axes = FALSE)
axis(1)
axis(2, at = (0:13)/13, labels = c("0", paste(1:13, 13, sep = "/")), las = 2 )
abline(h = cumsum(c(1,2,6,4)/13), lty = 3, col = "gray")
plot(e, add = TRUE, col = "red")

Answer 2

Your code generates some data d1 and associated numeric weights w1. These data are then treated as a reference population. The code takes a random sample d2 from this population of values d1, with sampling probabilities proportional to the associated weights w1. From the sample, you compute the weighted empirical distribution function f_d_2 of the sampled values d2, with weights inversely proportional to the original sampling probabilities. This function f_d_2 is a correct estimate of the original population distribution function, by the Horvitz-Thompson principle. But it's not exactly equal to the original population distribution, because it's a sample. The Kolmogorov-Smirnov test statistic should not be zero; it should be a small value.