PaperlayerView exact revision on arXiv

Source-linked AI summary

Determinantal point process models and statistical inference : Extended version

Frédéric Lavancier, Jesper Møller, Ege Rubak

arXiv:1205.4818v5math.ST

TL;DR

Statistical models and inference procedures for repulsive spatial point patterns remain underdeveloped, while Gibbs-process alternatives often have difficult likelihoods, moments, and simulation. The paper develops tractable parametric DPP models and inference tools, showing that they support closed-form moments and fast simulation, but remain less suitable for cases near Poisson or more repulsive than DPPs allow.

  • Problem

    Statistical models and inference procedures for repulsive spatial point patterns are largely unexplored, despite the need to describe datasets where nearby points repel.

  • Method

    The paper develops parametric stationary and inhomogeneous DPP models and uses their likelihood and moment properties for statistical inference, supported by spectral approximations and software.

  • Results

    DPP models provide closed-form moment expressions, tractable likelihood-based inference, and fast simulation for repulsive spatial point patterns.

  • Takeaways & Limitations

    DPPs provide parsimonious empirical models for comparing repulsive spatial point datasets through parameters, likelihoods, intensities, and pair correlation functions.

  • Takeaways & Limitations

    DPPs are less suitable near Poisson and cannot be as repulsive as Gibbs hard-core point processes.

Abstract

from arXiv · show

Statistical models and methods for determinantal point processes (DPPs) seem largely unexplored. We demonstrate that DPPs provide useful models for the description of spatial point pattern datasets where nearby points repel each other. Such data are usually modelled by Gibbs point processes, where the likelihood and moment expressions are intractable and simulations are time consuming. We exploit the appealing probabilistic properties of DPPs to develop parametric models, where the likelihood and moment expressions can be easily evaluated and realizations can be quickly simulated. We discuss how statistical inference is conducted using the likelihood or moment properties of DPP models, and we provide freely available software for simulation and statistical inference.

1 Introduction

The paper develops determinantal point processes (DPPs) as tractable statistical models for regular spatial point patterns, addressing a largely unexplored area of statistical inference. It establishes model, likelihood, moment, simulation, and inference tools while examining how much repulsiveness DPPs can represent.

  • DPPs model spatial point patterns in which nearby points repel, representing regularity rather than aggregation or clustering.
  • The paper surveys DPP definitions, existence conditions, moment properties, density expressions, and simulation procedures for statisticians.
  • It develops parametric model classes for stationary DPPs, later extending them to inhomogeneous DPPs, and studies their capacity to model repulsiveness.
  • The paper derives Fourier-series approximations for spectral representations to make simulation and inference feasible when explicit spectral decompositions are unavailable.
  • DPPs cannot match the repulsiveness of Gibbs hard-core processes, although the jinc-like DPP can fit quite regular datasets such as termite mounds.
  • DPPs offer tractable inference because moments are determinant-based, restricted processes have closed-form density expressions, and realizations can be simulated through mixtures of determinantal projection processes.

2 Definition, existence, simulation, and densities for determinantal point processes

DPPs define point-process product densities through determinants of a kernel, with Hermitian-kernel conditions ensuring repulsiveness and existence. Spectral representations support simulation and density evaluation on compact windows.

  • Definition: A DPP is defined by product densities equal to determinants of the kernel matrix evaluated at distinct points.The kernel must be non-negative definite; in the paper’s main setting it is a Hermitian complex covariance function.
  • Repulsiveness: For Hermitian kernels, the pair correlation satisfies g ≤ 1, and higher-order product densities decrease toward zero as points approach one another.This expresses the repulsive behavior of DPPs relative to a Poisson process.
  • Spectral representation: Mercer’s theorem provides a spectral representation on compact sets with unique eigenvalues and an orthonormal eigenfunction basis.Nonzero eigenvalues are real with finite multiplicity, and their only possible accumulation point is zero.
  • Existence: Under (C1), existence of DPP(C) is equivalent to condition (C2).The paper assumes both conditions for the remainder of its general DPP development.
  • Simulation: The implemented simulation algorithm is based on a mixture of determinantal projection point processes and generates the target compact-set DPP distribution.For a finite-rank kernel represented by n eigenfunctions, Algorithm 1 generates points distributed as the corresponding projection DPP.
  • Densities: When all eigenvalues are below one, the compact-set DPP has a density with respect to a unit-intensity homogeneous Poisson process.The density is expressed using an exponential factor and the determinant of the transformed kernel matrix.

3 Stationary models

The paper develops stationary, especially isotropic, DPP models through covariance and spectral constructions. These models expose a trade-off between intensity and correlation range while allowing flexible repulsiveness and practical approximations.

  • Stationarity and isotropy: Stationary DPP distributions are translation-invariant, and isotropic covariance functions make pair correlation depend only on interpoint distance.This permits use of pair-correlation and K-function procedures developed for spatial statistics.
  • Stationarity and isotropy: For isotropic models with decreasing covariance, pair correlation increases from zero toward one, motivating a distance-based range of correlation.The range r0 is a distance beyond which g0(r) is approximately one.
  • Spectral construction: The spectral approach characterizes admissible stationary models through a nonnegative spectral density ϕ with values bounded by one.The paper gives an equivalence between such spectral representations and the stationary DPP conditions under the stated integrability assumptions.
  • Model constraints: For fixed covariance-shape parameters, admissibility restricts intensity to 0 ≤ ρ ≤ ρmax, with ρmax decreasing as correlation range increases.This establishes a trade-off between how intense and how repulsive a stationary DPP can be.
  • Covariance models: The Whittle-Matérn family offers several pair-correlation shapes, while models with the same correlation range have maximal intensities of similar order.The latter observation indicates that interaction range strongly affects maximal permissible intensity.
  • Spectral approach: The power exponential spectral family contains the Gaussian model at ν = 2 and yields increasingly repulsive examples as ν increases toward infinity.At maximal permissible α, the paper approximates pair-correlation and L(r) − r functions because a closed form for the covariance is unavailable.

4 Approximations

The paper develops spectral approximations for DPP kernels, simulation, and densities, including periodic and border simulation methods for rectangular regions. The periodic approach is computationally efficient and empirically agrees closely with the border method, while truncated Fourier calculations support approximate density evaluation.

  • 4.1 Approximation of the kernel C: Spectral kernel approximations replace analytically unavailable representations using Fourier coefficients derived from the spectral density.For suitable covariance functions, the resulting approximate kernel is close to the true kernel, with accurate approximate pair correlation functions.
  • 4.2 Approximate simulation: The periodic method simulates an approximate DPP by imposing opposite-border interactions on S, producing low acceptance probabilities near boundaries.The method uses the approximation on S/2 to motivate simulation over the full region S.
  • 4.2 Approximate simulation: The periodic and border methods produce closely agreeing L(r) − r summaries across 1000 Gaussian-model realizations, suggesting nearly identical simulated DPPs.Figure 6 compares empirical means and pointwise quantiles for ρ = 100 and α = 0.05.
  • 4.2 Approximate simulation: The periodic method is preferred computationally: 1000 realizations were generated in approximately three minutes on a dual-core laptop.For circular covariance functions, the border method was exact and empirical summary-statistic plots showed almost no difference between methods.
  • 4.3 Approximation of the density: The density is approximated by f_app, with practical evaluation using truncated Fourier sums, symmetry reductions, or FFT-based interpolation.Likelihood inference based on f_app works well in the paper’s examples, whereas a convolution approximation can be poor in some situations.

5 Inference for stationary models

The paper develops likelihood- and moment-based inference for four stationary DPP model classes, using approximations that make estimation and model comparison feasible. Simulation studies and datasets show that sufficiently large truncations improve likelihood estimates, while Whittle–Matérn models often provide competitive fits with tractable moments.

  • Model setup: Stationary DPPs are modeled with Gaussian, Whittle–Matérn, Cauchy, and power exponential spectral parametric classes.The intensity is ρ, while θ parameterizes the correlation function.
  • Likelihood inference: Approximate maximum likelihood estimates θ by maximizing a truncated log-likelihood based on the determinant of an n × n kernel matrix.One-dimensional optimization can use a simple search; higher-dimensional optimization can use Nelder–Mead without explicit derivatives.
  • Likelihood approximation: N is increased until the approximate MLE stabilizes because it controls both spectral truncation and FFT grid resolution.The criterion S_N > 0.99ρ̂ alone may be insufficient for accurate likelihood approximation.
  • Intensity estimation: The largest relative difference between the non-parametric intensity estimate n/|S| and the MLE of ρ was 4% across real datasets.The authors therefore prefer ρ̂ = n/|S| for computational reasons, while noting that joint maximum likelihood gives similar estimates.
  • Simulation study: 500 simulated datasets showed that, when truncation is sufficiently large, MLE outperforms minimum contrast estimation because it has smaller biases and standard deviations.The likelihood approximation is less accurate for small N in the two Whittle–Matérn models because S_N converges more slowly.
  • Spanish towns dataset: In the Spanish towns analysis, the Whittle–Matérn model was preferred over comparable alternatives, and a simulation-based likelihood-ratio test rejected the Gaussian model with p-value 0.03.The fitted Whittle–Matérn model was judged to provide a good goodness-of-fit.
  • Model comparison: The Whittle–Matérn model can provide fewer parameters and direct access to intensity and pair correlation moments that require simulation for the Strauss hard-core model.The paper reports that it arguably provides a better fit than the Strauss hard-core model in the discussed example.
  • Termite mounds dataset: For the termite mounds dataset, the power exponential spectral model had the highest likelihood, while the Gaussian model was rejected with a simulation-based likelihood-ratio p-value of 2.6%.The fitted power exponential model was judged to provide a good fit despite lying on the border of the parameter space.

6 Inference for non-stationary models

The paper develops two strategies for fitting non-stationary DPPs: intensity-reweighted stationary models and transformations to stationary processes. Applications to Japanese pines and mucous membrane data illustrate model fitting, goodness-of-fit assessment, and practical limitations.

  • Model construction: Non-stationary DPPs can be constructed either by intensity-reweighting a stationary DPP or by transforming a stationary DPP.The first strategy is generally applicable, whereas the transformation approach is more dataset-specific.
  • 6.1 Second-order intensity-reweighted stationary models: For intensity-reweighted models, the intensity is parametrized, estimated by Poisson maximum likelihood, and used to reweight nonparametric summary statistics.A stationary isotropic dominating DPP provides the correlation structure, while the fitted intensity determines thinning probabilities.
  • 6.1 Second-order intensity-reweighted stationary models: For Japanese pines, log ρψ is modeled as a cubic polynomial in Cartesian coordinates, and the fitted Gaussian DPP provides an acceptable overall fit.The Gaussian model has α̂ = 0.226; intensity reestimation is required for each simulation when constructing valid envelopes.
  • 6.2 Inhomogeneity by transformation: A separable inhomogeneous intensity can be handled by transforming coordinates so the transformed process has unit intensity and can be fitted with stationary DPP models.The transformation maps the original process to a stationary-DPP fitting problem while preserving a DPP representation.
  • 6.2 Inhomogeneity by transformation: The transformed mucous membrane process is not second-order intensity-reweighted stationary, because its pair correlation depends on transformed coordinate differences.Its pair correlation is expressed as g(x,y) = 1 − |C_Y(T(y) − T(x))|^2.
  • 6.2 Inhomogeneity by transformation: For the mucous membrane data, the intensity is estimated piecewise on nine intervals, transformed, and then modeled using four stationary parametric DPP classes.The fitted power exponential model has the highest likelihood, but it is only slightly above the Gaussian model; the Gaussian is preferred despite a 5% lack-of-fit indication.

7 Concluding remarks

The paper presents DPPs as parsimonious empirical models for repulsive spatial patterns, with tractable likelihoods, moments, and simulation. It also identifies boundaries: DPPs cannot match hard-core repulsion, near-Poisson models are computationally difficult, and several extensions remain open.

  • Contributions: The paper introduces parametric DPP models and approximations that make likelihood evaluation and simulation practical.Its applications demonstrate likelihood-based and moment-based inference for spatial point patterns.
  • Interpretation: DPPs are intended as parsimonious empirical models for comparing datasets through fitted parameters, likelihoods, intensities, pair correlations, and other summaries.The authors explicitly distinguish this empirical role from modeling a physical mechanism generating the data.
  • Repulsiveness: At fixed intensity, stronger stationary-DPP repulsion is constrained, with the most repulsive planar model having a jinc-like kernel.The paper proposes μ as a rough measure of repulsiveness for stationary processes whose pair correlation does not exceed one.
  • Repulsiveness: DPPs cannot be as repulsive as Gibbs hard-core processes, although jinc-like DPPs can model highly regular datasets such as termite mounds.Compared with a fitted Strauss process, the jinc-like DPP offers closed-form moments and much faster simulation.
  • Computational limitations: Models close to Poisson are difficult to estimate and simulate because satisfactory spectral truncations require very large N.Weakly repulsive Gibbs models can therefore become competitive in this regime.
  • Open problems: The presented approximations are designed for rectangular observation windows, while fitting DPPs in non-rectangular windows remains unresolved.The paper identifies this as an area requiring further clarification.
  • Open problems: Future work includes DPPs on other spaces, multivariate and marked models, preferential sampling, incomplete observation, and space-time settings.The paper notes that inference for such extensions remains challenging.

Appendices

The appendices establish that smooth transformations and independent thinning preserve the DPP class. These closure properties support the paper’s constructions for changing domains and modeling inhomogeneity.

  • Transformations: A diffeomorphic transformation maps a DPP on one Borel set to a DPP on the transformed set.The transformed kernel is determined using the inverse transformation and its Jacobian.
  • Thinning: Independent thinning of a DPP with location-dependent retention probabilities produces another DPP.The resulting kernel incorporates the retention probabilities at both locations.

B Proof of Theorem 2.3

This appendix develops the technical conditions and distributions underlying DPP existence and simulation. It links local spectral properties to existence and gives recursion and inversion procedures for the number of points.

  • Existence conditions: DPP existence on compact sets is characterized through spectral conditions on the kernel restricted to those sets.Under continuity and non-negative definiteness, the required local trace-class properties follow from the kernel’s diagonal integral.
  • Palm distributions: Reduced Palm distributions of a DPP are themselves DPPs, with product densities related to higher-order product densities of the original process.The relationship is expressed by dividing the higher-order product density by the intensity.
  • Point-count distribution: The distribution of the number of points is computed from kernel eigenvalues using a recursion after the minimum feasible count m′.The probability mass below m′ is zero, and the remaining probabilities can be calculated recursively.
  • Simulation: The number of points can be simulated by inversion using the distribution function constructed from the recursively computed probabilities.A uniform random variable is mapped to the smallest count whose cumulative probability reaches it.

E Proof of Theorem 2.7 and related remarks

The appendix proves that the sequential sampling construction yields valid densities and produces an n-point DPP with joint density proportional to det[K]. It also explains numerically stable simulation and rejection-sampling considerations.

  • Proof of Theorem 2.7: The induction establishes that each p_i is a probability density and that the generated vectors remain linearly independent almost surely.The construction uses orthogonal projections onto complements of spans formed by previously generated vectors.
  • Proof of Theorem 2.7: The product of sequential projection norms equals the squared parallelepiped volume and therefore det[K](x_1, ..., x_n).This identity remains valid when the vectors are linearly dependent because both determinants vanish.
  • Proof of Theorem 2.7: The resulting n-point process has n-th order joint intensity ρ^(n)(x_1, ..., x_n) = det[K](x_1, ..., x_n).The fixed number of points converts the joint density into the n-th order product density.
  • Related remarks: Gram-Schmidt supplies orthonormal bases for the relevant subspaces, while projection norms provide an alternative interpretation of the sequential densities.For each i, i p_i(x) is the squared norm of the projection of v(x) onto H_i^⊥.
  • Related remarks: Recursive matrix multiplication can be numerically unstable, whereas Algorithm 1 calculates p_i(x) straightforwardly and stably.Uniform rejection sampling may become inefficient for small i, although the density computation remains fast in practice.
  • Related remarks: Explicit upper bounds on p_i yield easy-to-simulate instrumental densities for rejection sampling, including stepwise linear forms in one dimension.In two dimensions, the resulting instrumental density is a stepwise polynomial, with further upper-bounding potentially needed.

G Proof of Theorem 2.8 and related remarks

The appendix connects finite DPPs with conditional intensities, Poisson-process couplings, and repulsiveness. It also identifies restrictive conditions needed to extend these results globally.

  • Proof of Theorem 2.8: The n-th order product density can be represented as E f(Y ∪ {x_1, ..., x_n}), linking DPP densities to a Poisson-process functional.The cited result is attributed to prior work on DPP representations.
  • Related remarks: The hereditary density permits a Papangelou conditional intensity defined by a ratio of determinants involving ˜C.The ratio is taken as zero when both numerator and denominator are zero.
  • Proof of Theorem 2.8: The DPP’s monotonicity property confirms repulsiveness, while the associated conditional-intensity property establishes local stability.These properties are stated for the finite process X_S.
  • Related remarks: The finite DPP X_S can be coupled as a dependent thinning of a Poisson process with intensity function ˜C(u,u).The construction extends through increasing bounded windows to realize the full process as a dependent thinning of a Poisson process.
  • Related remarks: Global Papangelou intensities and reduced Palm distributions require finite-range conditions and sufficiently small C, which are particularly restrictive when d ≥2.The finite-process results themselves do not remove these scope restrictions.

H Proof of Proposition 3.1

The proof characterizes the spectrum of a stationary DPP’s integral operator through Fourier analysis. The operator spectrum equals the essential image of the spectral density multiplier.

  • Proof of Proposition 3.1: The eigenvalues and eigenfunctions in the kernel representation correspond to those of the integral operator T_S.This identifies the operator-theoretic objects used in the stationary analysis.
  • Proof of Proposition 3.1: For stationary kernels, extending functions by zero converts T_S into a convolution operator represented in Fourier space by multiplication with ϕ.The representation uses the Fourier transform and its inverse.
  • Proof of Proposition 3.1: The spectrum of T_S equals the spectrum of the multiplication operator Q_{ϕ,S}, namely the essential image of ϕ_S.This follows from the unitarity of the Fourier operator.

I Proof of Corollary 3.3

The proof establishes equivalence between spectral and kernel conditions for stationary DPPs using Fourier analysis, Bochner’s theorem, and operator bounds.

  • Proof of Corollary 3.3: A nonnegative integrable spectral density ϕ yields a continuous, positive-definite kernel C_0 through inverse Fourier transformation.Parseval’s identity also gives C_0 ∈ L2(R^d).
  • Proof of Corollary 3.3: Conversely, a continuous positive-definite C_0 in L2(R^d) has a nonnegative integrable Fourier density ϕ, bounded above by one.The nonnegativity follows from Bochner’s theorem, while the upper bound follows from Proposition 5.1.

J Quantifying and comparing repulsiveness

The paper compares DPP repulsiveness using K-functions and a rough scalar measure, then characterizes how model parameters control repulsion and approximation error.

  • Repulsiveness criteria: K-functions compare two stationary DPPs with equal intensity: smaller K(r) at every r indicates stronger repulsiveness.For isotropic pair-correlation functions, the criterion is consistent with comparing their corresponding pair-correlation behavior.
  • Parametric comparisons: At fixed ν, increasing α increases repulsiveness within the Gaussian, Whittle–Matérn, and Cauchy model classes, but lowers their maximal intensity.At α = αmax, Whittle–Matérn and Cauchy repulsiveness increases with ν and approaches the Gaussian case.
  • Parametric comparisons: The power exponential spectral model contains the Gaussian model at ν = 2 and approaches increasingly repulsive stationary DPPs as ν grows toward infinity.Its limiting case is the most repulsive stationary DPP, whose spectral density is an indicator function on a set of volume ρ.
  • Parametric comparisons: The K-function criterion cannot always order Whittle–Matérn and Cauchy models, motivating the supplementary scalar measure µ for broader comparisons.The measure is proposed for stationary point processes with constant intensity and translation-invariant pair correlation, when the integral exists.
  • Repulsiveness criteria: The scalar measure µ equals the limiting difference between expected point counts under the ordinary and reduced Palm distributions.For stationary processes, 0 ≤ µ ≤ 1; the Poisson process has µ = 0, while DPPs satisfy g ≤ 1 and therefore µ ≥ 0.
  • Approximation error: For the Whittle–Matérn kernel approximation, the error bound is small for reasonable ρ, ν, and α satisfying the model’s intensity condition.The bound becomes an equality when d = 1 and ν = 1/2; for higher dimensions, an upper bound is established using Bessel-function estimates.

L.2 Examples

The examples develop practical density approximations for likelihood computation, comparing convolution and periodic methods for Gaussian and related covariance models. Their convergence depends strongly on dimension and intensity, while convolution is faster but can be biased near the intensity boundary.

  • Approximation construction: The density approximation is implemented by truncating the infinite sums in the convolution and periodic expansions.The approximations require closed-form convolution powers for Gaussian and Whittle–Matérn models, but none is available for the Cauchy model.
  • Convergence: For Gaussian and Whittle–Matérn models, h⋆k(0) decays as k^-d/2, so convergence depends crucially on dimension d and intensity ρ.When d < 3, convergence requires ρ < ρmax and becomes slow as ρ approaches ρmax.
  • Gaussian example: Figure 17 compares convolution and periodic density approximations for the Gaussian log-likelihood as α varies, using a simulated unit-square dataset with ρ = 200 and α = 0.02.The dataset contains 213 observed points, and ρ is fixed at the observed intensity while α varies over (0, αmax), with αmax = 0.39.
  • Gaussian example: The periodic approximation gives effectively unbiased estimates, whereas the convolution approximation produces positively biased α estimates and often reaches αhat = αmax.This difference agrees with the likelihood comparison in Figure 17.
  • Gaussian example: When α < αmax/2, the two approximations are very similar and the convolution method converges rapidly.In this range, N = 10 is sufficient for stable results, compared with N = 512 for the periodic approximation, making convolution computationally faster.
Loading 1205.4818v5…