library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.univar
#> spatstat.univar 3.1-4
#> Loading required package: spatstat.geom
#> spatstat.geom 3.5-0.008
#> Loading required package: spatstat.random
#> spatstat.random 3.4-1
#> Loading required package: spatstat.explore
#> Loading required package: nlme
#> spatstat.explore 3.5-2
#> Loading required package: spatstat.model
#> Loading required package: rpart
#> spatstat.model 3.4-0.004
#> Loading required package: spatstat.linnet
#> spatstat.linnet 3.3-1
#>
#> spatstat 3.4-0.001
#> For an introduction to spatstat, type 'beginner'
# Define an intensity function in the unit square having the following expression:
# lambda(x, y) = 55 * (2x + 2y)
my_f <- \(x, y) 55 * (2 * x + 2 * y)
# Graphical representation
plot(as.im(my_f, W = owin()))
# Simulate points in the unit square according to the aforementioned intensity function
set.seed(1)
my_pp <- rpoispp(my_f, win = owin())
# Estimate the density of the simulated points using kernel density estimation
my_kd1 <- density(my_pp, sigma = 0.3)
plot(my_kd1)
points(my_pp, pch = 16)
# Simulate a point pattern from the estimated intensity
sim_kd1 <- rpoispp(my_kd1)
# Compute quadrat counts (4x4 grid) for the observed and simulated patterns
(qc_obs <- quadratcount(my_pp, nx = 4, ny = 4))
#> x
#> y [0,0.25] (0.25,0.5] (0.5,0.75] (0.75,1]
#> (0.75,1] 5 9 6 13
#> (0.5,0.75] 2 5 8 5
#> (0.25,0.5] 2 7 11 14
#> [0,0.25] 4 3 4 3
qc_sim <- quadratcount(sim_kd1, nx = 4, ny = 4)
# Compute a chi2 index
sum((qc_obs - qc_sim) ^ 2 / qc_sim)
#> [1] 49.14069
# Define another estimator using an extremely over-smoothing bandwidth
my_kd2 <- density(my_pp, sigma = 2)
plot(my_kd2)
points(my_pp, pch = 16)
# Compare the two estimates
solist(my_kd1, my_kd2) |> plot(equal.ribbon = TRUE, main = "")
# Repeat the previous process
sim_kd2 <- rpoispp(my_kd2)
qc_sim <- quadratcount(sim_kd2, nx = 4, ny = 4)
sum((qc_obs - qc_sim) ^ 2 / qc_sim)
#> [1] 57.07063
# The chi2 index is much larger for the over-smoothed estimate ---> we prefer the first one. Created on 2025-10-28 with reprex v2.1.1.9000