library(tidyverse)
library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.univar
#> spatstat.univar 3.1-1
#> Loading required package: spatstat.geom
#> spatstat.geom 3.3-5
#> Loading required package: spatstat.random
#> spatstat.random 3.3-2
#> Loading required package: spatstat.explore
#> Loading required package: nlme
#>
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#>
#> collapse
#> spatstat.explore 3.3-4
#> Loading required package: spatstat.model
#> Loading required package: rpart
#> spatstat.model 3.3-4
#> Loading required package: spatstat.linnet
#> spatstat.linnet 3.2-5
#>
#> spatstat 3.3-1
#> For an introduction to spatstat, type 'beginner'
set.seed(20250216)Thing 1 to notice:
when lambda is small enough, x can (w.p.a. 1) only obtain 0 and 1.
The probability you get a 1 is
poisson_small_lambda_is_approx_bernoulli <- function(lambda = 0.001) {
if (lambda < 0.00001) {
stop("lambda is too small for this numerical simulation to be precise")
}
B <- max(10000, 10 * (1 / lambda))
x <- rpois(B, lambda = lambda)
pct <- table(x, dnn = NULL) / B
cat(sprintf(
"Lambda: %s\nFreq. of draws: \n",
lambda
))
print(pct)
}
poisson_small_lambda_is_approx_bernoulli(lambda = 0.1)
#> Lambda: 0.1
#> Freq. of draws:
#> 0 1 2 3
#> 0.9054 0.0891 0.0050 0.0005
poisson_small_lambda_is_approx_bernoulli(lambda = 0.001)
#> Lambda: 0.001
#> Freq. of draws:
#> 0 1
#> 0.9986 0.0014
poisson_small_lambda_is_approx_bernoulli(lambda = 0.0001)
#> Lambda: 1e-04
#> Freq. of draws:
#> 0 1
#> 0.99991 0.00009
poisson_small_lambda_is_approx_bernoulli(lambda = 0.00001)
#> Lambda: 1e-05
#> Freq. of draws:
#> 0 1
#> 0.999991 0.000008Now, we will try to show that if a point process satisfies
- independent (across quadrats), and
- homogeneous intensity,
$\lambda$ , i.e. for all$B$ $\mathbb{E}(n(X \cap B)) = \lambda * \vert B \vert$
Then, it must be a Poisson point process:
To show this, I’m going to make a bunch of small cells
and put a point randomly in the cell if runif() is larger
than
Then, I will use big cells and count the number of points in each cell.
The counts should be distributed Poisson with intensity
This should be distributed Poisson with intensity
# Poisson point process intensity parameter
# Intensity per squared unit area
lambda <- 500
width <- 1
height <- 1
area <- width * height
mean_N <- lambda * area
create_cells <- function(lambda, width, height, nx, ny) {
total_area <- width * height
xmin <- seq(0, width, length.out = (nx + 1))[-(nx + 1)]
ymin <- seq(0, width, length.out = (ny + 1))[-(ny + 1)]
cells <- expand.grid(xmin = xmin, ymin = ymin)
cells$xmax <- cells$xmin + width / nx
cells$ymax <- cells$ymin + height / ny
cells$cell_area <- (cells$xmax - cells$xmin) * (cells$ymax - cells$ymin)
cells$cell_intensity <- lambda * cells$cell_area
return(cells)
}
# small boxes
cells_small <- create_cells(lambda, width, height, nx = 100, ny = 100)
# large boxes
nx <- 5
ny <- 5
draws <- c(replicate(10000, {
# Poisson is approximately bernoulli with probability `cell_intensity`
cells_small$has_point <-
runif(nrow(cells_small)) <= cells_small$cell_intensity
# for cells with points, grab a random location on the cell
points <- cells_small[cells_small$has_point == TRUE, ]
points$y <- points$ymin + runif(nrow(points)) * (points$ymax - points$ymin)
points$x <- points$xmin + runif(nrow(points)) * (points$xmax - points$xmin)
# Count cells in cells_large
large_cell_counts <- quadratcount(
ppp(
points$x, points$y,
window = owin(xrange = c(0, width), yrange = c(0, height))
),
ny = 5, nx = 5
)
large_cell_counts <- c(large_cell_counts)
return(large_cell_counts)
}))
# Compare with poisson distribution with area
draws_poisson <- rpois(10000, lambda = lambda * area / (nx * ny))library(tinyplot)
tinytheme("clean2")
df <- rbind(
data.frame(type = "Large Cell Counts", draws = c(draws)),
data.frame(type = "Poisson", draws = draws_poisson)
)
tinyplot(
~draws,
facet = ~type,
data = df,
type = "histogram",
freq = FALSE,
breaks = -0.5 + 0:(max(df$draws) + 1),
xlab = ""
)