Edit on GitHub

Spatial point processes

Partial rejection sampling (PRS) can be applied to generate exact samples from pairwise Gibbs point processes with finite interaction range.

For more details on spatial point processes please refer to Jesper. Møller , Rasmus Plenge. Waagepetersen (2004) and references therein.

Poisson point process

PartialRejectionSampling.HomogeneousPoissonPointProcess — Type
HomogeneousPoissonPointProcess{T<:Vector{Float64}} <: AbstractSpatialPointProcess{T}

Homegeneous Poisson point process with intensity $\beta > 0$, denoted $\operatorname{Poisson}(\beta)$.

$\operatorname{Poisson}(\beta)$ has density (w.r.t. the homogenous Poisson point process with unit intensity $\operatorname{Poisson}(1)$) proportional to

\[ \prod_{x \in X} \beta.\]

source
PartialRejectionSampling.HomogeneousPoissonPointProcess — Type
HomogeneousPoissonPointProcess(
    β::Real,
    window::Union{Nothin,AbstractSpatialWindow}=nothing
)

Construct a PRS.HomogeneousPoissonPointProcess with intensity β restricted to window.

Default window (window=nothing) is PRS.SquareWindow().

using PartialRejectionSampling

β = 40
win = PRS.SquareWindow(zeros(2), 1)
PRS.HomogeneousPoissonPointProcess(β, win)

# output

HomogeneousPoissonPointProcess{Array{Float64,1}}
- β = 40.0
- window = SquareWindow [0.0, 1.0]^2
source
PartialRejectionSampling.generate_sample — Method
generate_sample(
    [rng::Random.AbstractRNG,]
    pp::HomogeneousPoissonPointProcess{Vector{T}},
    win::Union{Nothing,AbstractWindow}=nothing
)::Matrix{T} where {T<:Float64}

Generate an exact sample from PRS.HomogeneousPoissonPointProcess on window win. Sampled points are stored as columns of the output matrix.

Default window (win=nothing) is window(pp).

source
PartialRejectionSampling.generate_sample_poisson_union_balls — Function
generate_sample_poisson_union_balls(
    [rng::Random.AbstractRNG,]
    β::Real,
    centers::Matrix,
    radius::Real,
    win::Union{Nothing,AbstractWindow}=nothing
)::Matrix{Float64}

Generate an exact sample from a homogenous PRS.HomogeneousPoissonPointProcess Poisson(β) on $\bigcup_{i} B(c_i, r)$ (union of balls centered at $c_i$ with the same radius $r$).

If win ≂̸ nothing, returns the points falling in win.

Hint

Use the independence property of the Poisson point process on disjoint subsets in order to

  • Sample from Poisson(β) on $B(c_1, r)$,
  • Sample from Poisson(β) on $B(c_2, r) \setminus B(c_1, r)$,
  • Sample from Poisson(β) on $B(c_j, r) \setminus \bigcup_{i<j} B(c_i, r)$,
  • ...
source

Hard core point process

PartialRejectionSampling.HardCorePointProcess — Type
HardCorePointProcess{T<:Vector{Float64}} <: AbstractSpatialPointProcess{T}

Spatial point process with density (w.r.t. the homogenous Poisson point process with unit intensity) proportional to

\[ \prod_{x \in X} \beta \prod_{\{x, y\} \subseteq X} 1_{ \left\| x - y \right\|_2 > r },\]

where $\beta > 0$ is called the background intensity and $r > 0$ the interaction range.

PRS.HardCorePointProcess can be viewed as

See also

Example

A realization for $\beta=38$ and $r=0.1$ on $[0, 1]^2$, where points are marked with a circle of radius $r/2$.

assets/hard_core_spatial.png

source
PartialRejectionSampling.HardCorePointProcess — Type
HardCorePointProcess(
    β::Real,
    r::Real,
    window::Union{Nothing,AbstractSpatialWindow}=nothing
)

Construct a PRS.HardCorePointProcess with intensity β and interaction range r, restricted to window.

Default window (window=nothing) is PRS.SquareWindow().

using PartialRejectionSampling

β, r = 40, 0.05
win = PRS.SquareWindow(zeros(2), 1)
PRS.HardCorePointProcess(β, r, win)

# output

HardCorePointProcess{Array{Float64,1}}
- β = 40.0
- r = 0.05
- window = SquareWindow [0.0, 1.0]^2
source
PartialRejectionSampling.generate_sample_prs — Method
generate_sample_prs(
    [rng::Random.AbstractRNG,]
    pp::HardCorePointProcess{T},
    win::Union{Nothing,AbstractWindow}=nothing
)::Vector{T} where {T}

Sample from PRS.HardCorePointProcess using Partial Rejection Sampling (PRS) of Heng Guo , Mark Jerrum (2018).

Default window (win=nothing) is window(pp)=pp.window.

See also

Example

A illustration of the procedure for $\beta=38$ and $r=0.1$ on $[0, 1]^2$ where points are marked with a circle of radius $r/2$. Red disks highlight bad regions and orange disk, with radius $r$, describe the resampling region.

assets/hard_core_spatial_prs.gif

source
PartialRejectionSampling.gibbs_interaction — Method
gibbs_interaction(
    pp::HardCorePointProcess{T},
    cell1::SpatialCellGridPRS{T},
    cell2::SpatialCellGridPRS{T}
)::Real where {T}

Compute the pairwise Gibbs interaction for a PRS.HardCorePointProcess between $C_1$=cell1 and $C_2$=cell2,

\[ \prod_{(x, y) \in C_1 \times C_2} 1_{ \left\|x - y\right\|_2 \leq r }.\]

source

Strauss point process

PartialRejectionSampling.StraussPointProcess — Type
StraussPointProcess{T<:Vector{Float64}} <: AbstractSpatialPointProcess{T}

Spatial point process with density (w.r.t. the homogenous Poisson point process with unit intensity) proportional to

\[ \prod_{x \in X} \beta \prod_{\{x, y\} \subseteq X} \gamma^{ 1_{ \left\| x - y \right\|_2 \leq r } } = \beta^{|X|} \gamma^{|\{\{x, y\} \subseteq X ~;~ \left\| x - y \right\|_2 \leq r\}|},\]

with intensity $\beta > 0$, interaction coefficient $0\leq \gamma\leq 1$ and interaction range $r > 0$.

See also

source
PartialRejectionSampling.StraussPointProcess — Type
StraussPointProcess(
    β::Real,
    γ::Real,
    r::Real,
    window::Union{Nothing,AbstractSpatialWindow}=nothing
)

Construct a PRS.StraussPointProcess with intensity β, interaction coefficient γ, and interaction range r, restricted to window.

Default window (window=nothing) is PRS.SquareWindow().

using PartialRejectionSampling

β, γ, r = 2.0, 0.2, 0.7
win = PRS.SquareWindow([0.0, 0.0], 10.0)
PRS.StraussPointProcess(β, γ, r, win)

# output

StraussPointProcess{Array{Float64,1}}
- β = 2.0
- γ = 0.2
- r = 0.7
- window = SquareWindow [0.0, 10.0]^2

Example

A illustration of the procedure for $\beta=78, \gamma=0.1$ and $r=0.07$ on $[0, 1]^2$ where points are marked with a circle of radius $r/2$.

source
PartialRejectionSampling.gibbs_interaction — Method
gibbs_interaction(
    pp::StraussPointProcess{T},
    cell1::SpatialCellGridPRS{T},
    cell2::SpatialCellGridPRS{T}
)::Real where {T}

Compute the pairwise Gibbs interaction for a PRS.StraussPointProcess between $C_1$=cell1 and $C_2$=cell2,

\[ \prod_{(x, y) \in C_1 \times C_2} \gamma^{ 1_{ \left\|x - y\right\|_2 \leq r } }.\]

source

Windows

PartialRejectionSampling.RectangleWindow — Type
RectangleWindow{T<:Float64} <: AbstractRectangleWindow{T}

Structure representing a hyperrectangle $\prod_i [c_i, c_i + w_i]$, with fields

  • c lower left corner of the hyperrectangle
  • w width vector of the hyperrectangle along each coordinate
source
Base.in — Method
Base.in(
    x::AbstractVector,
    win::AbstractRectangleWindow
)

Check if $x \in B(c, r)$, i.e., $\left\| x - c \right\| \leq r$.

source
Base.in — Method
Base.in(
    x::AbstractVector,
    win::AbstractRectangleWindow
)

Check if $x \in \prod_{i} [c_i, c_i + w_i]$.

source
Base.rand — Method
Base.rand(
    [rng::Random.AbstractRNG,]
    win::BallWindow{T},
    n::Int
)::Matrix{T} where {T}

Sample n points uniformly at random in window win.

source
Base.rand — Method
Base.rand(
    [rng::Random.AbstractRNG,]
    win::BallWindow{T}
)::Vector{T} where {T}

Sample uniformly at random in window win.

source
Base.rand — Method
Base.rand(
    [rng::Random.AbstractRNG,]
    win::AbstractRectangleWindow{T},
    n::Int
)::Matrix{T} where {T}

Sample n points uniformly at random in window win.

source
Base.rand — Method
Base.rand(
    [rng::Random.AbstractRNG,]
    win::AbstractRectangleWindow{T}
)::Vector{T} where {T}

Sample uniformly at random in window win.

source