Now I have a clear picture. Here's the full mathematical framing of the simplest LGCP (Log-Gaussian Cox Process) model. --- ## The Simplest LGCP: Math and Framing ### The Two-Layer Story A Log-Gaussian Cox Process is a **doubly stochastic** point process. There are two layers of randomness: 1. **Layer 1 (Latent):** The log-intensity is a Gaussian Process (random function) 2. **Layer 2 (Observable):** Conditional on that intensity, points are generated as a Poisson process --- ### Formal Definition Let \( X \subseteq \mathbb{R}^d \) be the domain (e.g. an interval in 1D, or a region in 2D). **Step 1: The intensity field** Define a latent Gaussian process over the domain: \[ f(s) \sim \mathcal{GP}\big(\mu(s),\, K(s, s')\big), \qquad s \in X \] where: - \(\mu(s)\) is the **mean function** (often set to a constant \(\mu\)) - \(K(s,s')\) is the **covariance kernel** (e.g. the squared exponential / RBF kernel) To ensure the intensity \(\lambda(s)\) is strictly positive, exponentiate: \[ \lambda(s) = \exp\big(f(s)\big) \] So the **log-intensity is a Gaussian process** — hence the name "Log-Gaussian Cox Process." **Step 2: The point process** Conditional on a realization of \(\lambda(s)\), the points \(\{s_1, \dots, s_N\}\) are drawn from an **inhomogeneous Poisson process** with intensity \(\lambda(s)\). This means: - For any subset \(A \subseteq X\), the number of points \(Y_A\) falling in \(A\) follows: \[ Y_A \mid \lambda(\cdot) \sim \text{Poisson}\!\left( \int_A \lambda(s)\, ds \right) \] - The total number of points \(N\) in the whole domain is: \[ N \mid \lambda(\cdot) \sim \text{Poisson}\!\left( \int_X \lambda(s)\, ds \right) \] - Given \(N\), the points are independent and identically distributed with density proportional to \(\lambda(s)\): \[ p(s_1, \dots, s_N \mid N, \lambda(\cdot)) = \prod_{i=1}^N \frac{\lambda(s_i)}{\int_X \lambda(s)\, ds} \] --- ### The Simplest Special Case (1D, Constant Mean, RBF Kernel) This is the "hello world" of LGCPs. On a 1D domain (say \(t \in [0, T]\)): \[ \begin{aligned} \mu(t) &= \mu_0 \quad (\text{constant}) \\[4pt] K(t, t') &= \sigma^2 \exp\!\left(-\frac{(t-t')^2}{2\ell^2}\right) \quad (\text{RBF kernel}) \\[4pt] f(t) &\sim \mathcal{GP}(\mu_0, K) \\[4pt] \lambda(t) &= \exp(f(t)) \\[4pt] N \mid \lambda(\cdot) &\sim \text{Poisson}\!\left( \int_0^T \lambda(t)\, dt \right) \end{aligned} \] **Hyperparameters:** - \(\mu_0\) — baseline log-intensity (controls average number of points) - \(\sigma^2\) — marginal variance of the GP (controls how much the intensity fluctuates) - \(\ell\) — lengthscale (controls how smoothly the intensity varies in space/time) --- ### How to Simulate (Thinning Algorithm) A simple recipe from the blog post: 1. **Sample the GP:** at \(n\) grid points \(\{x_1,\dots,x_n\}\), draw \(f \sim \mathcal{N}(0, K)\), set \(\lambda_i = \exp(f_i)\) 2. **Find the maximum intensity:** \(M = \max_i \lambda_i\) 3. **Generate candidate events:** Draw \(N_{\text{cand}} \sim \text{Poisson}(M \cdot T)\) and place them uniformly 4. **Thin:** Accept each candidate at position \(t_i\) with probability \(\lambda(t_i) / M\) This yields a valid sample from the LGCP. --- ### Grid-Based Inference (Practical Fitting) In practice, since the GP is infinite-dimensional, we discretize the domain into \(M\) cells \(A_1,\dots,A_M\): \[ Y_j \mid \lambda_j \sim \text{Poisson}(\lambda_j \cdot |A_j|), \qquad \log\boldsymbol{\lambda} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{K}) \] where: - \(Y_j\) is the observed count in cell \(j\) - \(|A_j|\) is the area of cell \(j\) - \(\log\boldsymbol{\lambda} = [\log\lambda_1, \dots, \log\lambda_M]^\top\) - \(\mathbf{K}_{ij} = K(s_i, s_j)\) where \(s_i\) is the centroid of cell \(i\) This is a **Poisson GLMM with a GP random effect** — and can be fit via MCMC (PyMC, Stan), INLA, or variational inference. --- ### Summary Table | Component | Mathematical object | Role | |---|---|---| | **GP prior** | \(f(s) \sim \mathcal{GP}(\mu, K)\) | Models latent log-intensity as a random smooth function | | **Exponential link** | \(\lambda(s) = \exp(f(s))\) | Ensures positive intensity | | **Poisson likelihood** | \(Y_A \mid \lambda \sim \text{Poi}(\int_A \lambda)\) | Generates observed point counts | | **Kernel (RBF)** | \(K(s,s') = \sigma^2 \exp(-\|s-s'\|^2/2\ell^2)\) | Controls smoothness and correlation length of intensity | That's the simplest complete formulation of a Cox process modeled as a Gaussian process — the **Log-Gaussian Cox Process (LGCP)**.