Edit on GitHub

Sampling methods

Below is the list of point processes which can be sampled exactly using PartialRejectionSampling.jl

Default exact sampler

PartialRejectionSampling.generate_sample — Function
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
generate_sample(
    [rng::Random.AbstractRNG,]
    pp::HardCorePointProcess{T},
    win::Union{Nothing,AbstractWindow}=nothing
)::Vector{T} where {T}

Generate an exact sample from the PRS.HardCorePointProcess.

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

Default sampler is PRS.generate_sample_prs.

See also

source
generate_sample(
    [rng::Random.AbstractRNG,]
    pp::StraussPointProcess{T},
    win::Union{Nothing,AbstractWindow}=nothing
)::Vector{T} where {T}

Genererate an exact sample from PRS.StraussPointProcess.

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

Default sampler is PRS.generate_sample_dcftp.

See also

source
generate_sample(
    [rng::Random.AbstractRNG,]
    pp::Ising{T},
    idx::T
)::T where {T<:Int}

Generate an exact sample from the marginal distribution of state $i=$ idx of pp.

More specifically,

\[ x_i \sim \operatorname{Bernoulli}_{-1, 1} (\sigma(h_i)),\]

where $\sigma$ denotes the sigmoid function.

source
generate_sample(
    [rng::Random.AbstractRNG,]
    pp::Ising;
    win::GraphNode
)

Generate an exact sample from the marginal distribution of pp at state indexed by win.idx.

source
generate_sample(
    [rng=Random.AbstractRNG,]
    pp::HardCoreGraph{T}
)::Vector{T} where {T}

Generate an exact sample from the PRS.SinkFreeGraph.

Default sampler is PRS.generate_sample_prs.

source
generate_sample(
    [rng::Random.AbstractRNG,]
    pp::RootedSpanningForest
)

Generate an exact sample from the PRS.RootedSpanningForest.

Default sampler is PRS.generate_sample_prs.

source
generate_sample([rng::Random.AbstractRNG,] pp::SinkFreeGraph)

Generate an exact sample from the PRS.SinkFreeGraph.

Default sampler is PRS.generate_sample_prs.

source
generate_sample(
    [rng=Random.AbstractRNG,]
    pp::PatternFreeString{T},
    size::Int
)::T where {T<:String}

Generate a string uniformly at random among all strings made of characters from pp.alphabet with no occurence of pp.pattern.

Default sampler is PRS.generate_sample_prs.

source

Partial Rejection Sampling (PRS)

Heng Guo , Mark Jerrum , Jingcheng Liu (2019) developed the PRS methodology to generate exact samples from product distributions of the form $\otimes_{n=1}^{N} \mu_n$ subject to some constraints. It requires access to an exact sampler for each $\mu_n$ and an oracle to check the violation of the constraints.

Given an initial sample from $\otimes_{n=1}^{N} \mu_n$:

  • Vanilla rejection sampling resample all variables if any constraint is violated; until all constraints are satisfied,
  • Partial rejection sampling instead constructs a subset of variables to be resampled, starting from variables involved in violated constraints, and preserves the state of the variables outside of this resampling set; until all constraints are satisfied.

In both cases, the output sample is guaranteed to have the right distribution, i.e., the product distribution subject to the prescribed constraints.

PartialRejectionSampling.generate_sample_prs — Function
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
generate_sample_prs(
    [rng::Random.AbstractRNG,]
    pp::StraussPointProcess{T},
    win::Union{Nothing,AbstractWindow}=nothing
)::Vector{T} where {T}

Genererate an exact sample from PRS.StraussPointProcess using Partial Rejection Sampling (PRS).

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

Default sampler is PRS.generate_sample_grid_prs.

source
generate_sample_prs(
    [rng::Random.AbstractRNG,]
    pp::Ising
)

Generate an exact sample form pp using Partial Rejection Sampling (PRS), see Section 4.2 of Heng Guo , Mark Jerrum , Jingcheng Liu (2019).

Default sampler is PRS.generate_sample_grid_prs.

See also

source
generate_sample_prs(
    [rng=Random.AbstractRNG,]
    pp::HardCoreGraph{T}
)::Vector{T} where {T}

Sample from PRS.HardCoreGraph using Partial Rejection Sampling (PRS), see Section 7.2 of Heng Guo , Mark Jerrum , Jingcheng Liu (2019)

See also

Example

An illustration of the procedure on a $5\times 5$ grid graph.

  • red: variables involved in constraints violation (gray neighboring nodes)
  • orange: variables to be resampled (red nodes and their neighborhood)

assets/hard_core_graph.gif

source
generate_sample_prs(
    [rng::Random.AbstractRNG,]
    pp::RootedSpanningForest{T}
)::T where {T<:LG.SimpleDiGraph{Int64}}

Generate a rooted spanning forest of pp.graphwith prescribedpp.roots, uniformly at random among all rooted spanning forests rooted atpp.roots`, using Partial Rejection Sampling (PRS).

See also

Example

An illustration of the procedure on a $5\times 5$ grid graph, with roots=[13]

  • red: variables involved in constraints violation (edges forming a cycle)
  • orange: variables resampled (edges originating from a red/orange node)

assets/rooted_spanning_tree_prs.gif

source
generate_sample_prs(
    [rng::Random.AbstractRNG,]
    pp::SinkFreeGraph{T}
)::T where {T}

Generate an orientated version of pp.graph uniformly at random among all possible orientations conditioned on the absence of sinks, using Partial Rejection Sampling (PRS).

See also

Example

A illustration of the procedure on a $5 \times 5$ grid grah.

  • red: variables involved in constraints violation (orientation of the edges pointing to a sink node)
  • orange: variables resampled (orientation of the previous edges)

assets/sink_free_graph_prs.gif

source
generate_sample_prs(
    [rng::Random.AbstractRNG,]
    pp::PatternFreeString{T},
    size::Int
)::T where {T<:String}

Generate a string uniformly at random among all strings made of characters from pp.alphabet with no occurence of pp.pattern, using a tailored version of Partial Rejection Sampling (PRS).

using PartialRejectionSampling
pp = PRS.PatternFreeString("ATGTA", ["A", "C", "G", "T"])
PRS.generate_sample_prs(pp, 20)

See also

source

Grid Partial Rejection Sampling (grid PRS)

Sarat B. Moka , Dirk P. Kroese (2020) adapted the idea of Partial Rejection Sampling (PRS) (originally derived in the finite setting) to generate exact samples from Spatial point processes with finite range of interaction, noted $r > 0$.

The name of the method grid PRS results from the combination of the Partial Rejection Sampling (PRS) methodology with a partitioning of the domain of $\mathbb{R}^2$ where the target PRS.AbstractSpatialPointProcess – is divided into $N$ cells of type PRS.SpatialCellGridPRS with size $r \times r$,

To draw the correspondence with framework of PRS developed by Heng Guo , Mark Jerrum , Jingcheng Liu (2019), $\mu_n$ represents the target PRS.AbstractSpatialPointProcess restricted to the $n$-th cell, i.e., variables are point configurations.

Note that the grid PRS methodology requires an efficient sampling algorithm to generate exact samples on each cell involved in the partitioning of original domain, see, e.g., PRS.generate_sample_dcftp and PRS.generate_sample_prs.

See also closely related variants of grid PRS for

Base.iterate — Function
Base.iterate(cell::AbstractCellGridPRS, state=1) = iterate(cell.value, state)
source
PartialRejectionSampling.find_bad_cells_indices! — Method
find_bad_cells_indices!(
    [rng::Random.AbstractRNG,]
    g::SWG.SimpleWeightedGraph{T,U},
    cells::Vector{V},
    pp::AbstractPointProcess
)::Set{T} where {T,U,V<:AbstractCellGridPRS}

Identify bad events and return the corresponding cells' index. An event $\{i,j\}$ is said to be "bad"

\[ \left\{U_{ij} > \exp \left[ -\sum_{x \in C_i} \sum_{y \in C_j} V(x,y) \right] \right\}\]

where $U_{ij}$ is the weight of the edge $\{i,j\}$ in the interaction graph $g$ created by PRS.weighted_interaction_graph and $V$ the Gibbs potential discribing the pairwise Gibbs interaction of pp

Note when a bad event occurs, the corresponding $U_{ij}$ is resampled hence the "!"

This is a subroutine of PRS.generate_sample_grid_prs.

source
PartialRejectionSampling.find_cells_to_resample_indices! — Method
find_cells_to_resample_indices!(
    rng::Random.AbstractRNG,
    g::SWG.SimpleWeightedGraph{T,U},
    cells::Vector{V},
    pp::AbstractPointProcess
)::Set{T} where {T,U,V<:AbstractCellGridPRS}

Identify the set of events to be resampled as constructed by Algorithm 5 in Heng Guo , Mark Jerrum , Jingcheng Liu (2019) as part of the Partial Rejection Sampling (PRS) method. Return the indices of the variables (here cells) involved in the corresponding events.

This function is used as a subroutine of the grid PRS methodology of Sarat B. Moka , Dirk P. Kroese (2020), see PRS.generate_sample_grid_prs.

Note if the event associated to the edge $\{i,j\}$ of g is selected to be resampled, the uniform random variable encoded as the weight of the correspond edge is resampled (hence the "!")

This is a subroutine of PRS.generate_sample_grid_prs.

source
PartialRejectionSampling.is_inner_interaction_possible — Method
is_inner_interaction_possible(
    spp::AbstractSpatialPointProcess,
    cell_i::SpatialCellGridPRS,
    cell_j::SpatialCellGridPRS
) = false

Assume cell_i and cell_j are neighboring cells in the weighted interaction graph constructed by PRS.weighted_interaction_graph from spp and already identified in the set of variables to be resampled in PRS.generate_sample_grid_prs Since the configuration of points in the corresponding cells and the uniform random variable associated to the event \{i,j\} are considered fixed, there is no degree of freedom to make the interaction between cell_i and cell_j possible.

This is a subroutine of PRS.generate_sample_grid_prs.

source
PartialRejectionSampling.is_inner_interaction_possible — Method
is_inner_interaction_possible(
    ising::Ising{T},
    xᵢ::GraphCellGridPRS{T},
    xⱼ::GraphCellGridPRS{T}
) = true

Assume xᵢ and xⱼ are neighboring sites in the weighted interaction graph constructed by PRS.weighted_interaction_graph from ising.graph and already identified in the set of variables to be resampled PRS.generate_sample_grid_prs. Given the states of xᵢ and xⱼ, check whether a new assigment of $U_{ij}$ (the weight of edge $\{i,j\}$) can induce the bad event

\[ U_{ij} > \exp(J x_i x_j - |J|))\]

This is a subroutine of PRS.generate_sample_grid_prs.

source
PartialRejectionSampling.is_outer_interaction_possible — Method
is_outer_interaction_possible(
    spp::AbstractSpatialPointProcess,
    cell_i::SpatialCellGridPRS,
    cell_j::SpatialCellGridPRS
)::Bool

Assume cell_i and cell_j are neighboring cells in the weighted interaction graph constructed by PRS.weighted_interaction_graph from spp. Given the configuration of points in cell_i, check whether a realization of spp in cell_j can induce a bad event.

This is a subroutine of PRS.generate_sample_grid_prs.

source
PartialRejectionSampling.is_outer_interaction_possible — Method
is_outer_interaction_possible(
    ising::Ising{T},
    xᵢ::GraphCellGridPRS{T},
    xⱼ::GraphCellGridPRS{T}
) = true

Assume xᵢ and xⱼ are neighboring sites in the weighted interaction graph constructed by PRS.weighted_interaction_graph from ising.graph. Given the state of xᵢ, one can always find a new assigment of xⱼ and/or $U_{ij}$ (the weight of edge $\{i,j\}$) can induce the bad event

\[ U_{ij} > \exp(J x_i x_j - |J|))\]

This is a subroutine of PRS.generate_sample_grid_prs.

source
PartialRejectionSampling.weighted_interaction_graph — Method
weighted_interaction_graph(
    rng::Random.AbstractRNG,
    pp::AbstractSpatialPointProcess;
)::SWG.SimpleWeightedGraph

Construct the weighted interaction graph (King graph) used in PRS.generate_sample_grid_prs, to generate exact samples from PRS.AbstractSpatialPointProcess

The pp.window is divided into cells of length the interaction range pp.r. Each cell represents a vertex of the interaction (king) graph and each edge carries a uniform random varialble.

See also

source

Uniform sampling in spatial windows

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

Sample uniformly at random in window win.

source
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(
    [rng::Random.AbstractRNG,]
    win::BallWindow{T}
)::Vector{T} where {T}

Sample uniformly at random in window win.

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

Sample n points uniformly at random in window win.

source

Dominated Coupling From The Past (dCFTP)

Implementation of dominated Coupling From The Past (dCFTP) developed by Wilfrid Stephen Kendall , Jesper Møller (1999) and Wilfrid Stephen Kendall , Jesper Møller (2000)

See also

PartialRejectionSampling.backward_extend! — Method
backward_extend!(
    rng::Random.AbstractRNG,
    D::Set{T},          # Dominating process
    M::Vector{Float64}, # Marking process
    R::Vector{T},       # Recording process
    steps::StepRange,   # Number of backward steps
    birth_rate::Real,
    win::AbstractSpatialWindow
) where {T}

Sample from the dominating birth-death process backwards in time according to steps. The marks M and the points R which were added (uniform mark) / deleted (mark=0) along the run are recorded (pushfirst!). The final state of the dominating process is updated to D.

See also

source
PartialRejectionSampling.forward_coupling — Method
forward_coupling(
    D::Set{T},           # Dominating process
    M::Vector{Float64},  # Marking process
    R::Vector{T},        # Recording process
    pp::AbstractSpatialPointProcess{T},
    β::Real              # Upper bound on papangelou conditional intensity
) where {T}

Check if coalescence occured between the lower and upper bounding processes and return the state of the lower bounding process at time 0.

See also

source
PartialRejectionSampling.generate_sample_dcftp — Method
generate_sample_dcftp(
    [rng::Random.AbstractRNG,]
    pp::AbstractSpatialPointProcess{T};
    win::Union{Nothing,AbstractSpatialWindow}=nothing,
    n₀::Int=1
)::Vector{T} where {T}

Generate an exact sample from a spatial point process pp on window win using dominated coupling from the past.

  • Default window (win=nothing) is window(pp)=pp.window
  • Initial coalescence check performed after n₀ steps.
  • Seed or random number generator is passed via rng.

See also

source