And here is the R-code I came up with
`library(deldir) library(tidyverse) library(fields) #MPLE # irregular parameter k k <- 0.4 ## Generate dummy points 50X50. "RA" and "DE" are x and y coordinates dum.x <- seq(ramin, ramax, length = 50) dum.y <- seq(demin, demax, length = 50) dum <- expand.grid(dum.x, dum.y) colnames(dum) <- c("RA", "DE") ## Combine with data and specify which is data point and which is dummy, X is the point pattern to be fitted bind.x <- bind_rows(X, dum) %>% mutate(Ind = c(rep(1, nrow(X)), rep(0, nrow(dum)))) ## Calculate Quadrature weights using Voronoi cell area w <- deldir(bind.x$RA, bind.x$DE)$summary$dir.area ## Response y <- bind.x$Ind/w # the sum of distances between all pairs of points (the sufficient statistics) tmp <- cbind(bind.x$RA, bind.x$DE) t1 <- rdist(tmp)^(-2/k) t1[t1 == Inf] <- 0 t1 <- rowSums(t1) t <- -t1 # fit the model using quasipoisson regression fit <- glm(y ~ t, family = quasipoisson, weights = w) `
However, the fitted parameter for
t is negative which is obviously not a correct value for a softcore point process. Also, my point pattern is actually simulated from a softcore process so it does not make sense that the fitted parameter is negative. I tried my best to find any bugs in the code but I can't seem to find it. The only potential issue I see is that my sufficient statistics is extremely large (on the order of 10^14) which I fear may cause numerical issues. But the statistics are large because my observation window spans a very small unit and the average distance between a pair of points is around 0.006. So sufficient statistics based on this will certainly be very large and my intuition tells me that it should not cause a numerical problem and make the fitted parameter to be negative.
Can anybody help and check if my code is correct? Thanks very much!