Optimality of Poisson processes intensity learning with Gaussian processes

In this paper we provide theoretical support for the so-called"Sigmoidal Gaussian Cox Process"approach to learning the intensity of an inhomogeneous Poisson process on a $d$-dimensional domain. This method was proposed by Adams, Murray and MacKay (ICML, 2009), who developed a tractable computational approach and showed in simulation and real data experiments that it can work quite satisfactorily. The results presented in the present paper provide theoretical underpinning of the method. In particular, we show how to tune the priors on the hyper parameters of the model in order for the procedure to automatically adapt to the degree of smoothness of the unknown intensity and to achieve optimal convergence rates.


Introduction
Inhomogeneous Poisson processes are widely used models for count and point data in a variety of applied areas. A typical task in applications is to learn the underlying intensity of a Poisson process from a realised point pattern.
In this paper we consider nonparametric Bayesian approaches to this problem. These do not assume a specific parametric form of the intensity function and produce posterior distributions which do not only give an estimate of the intensity, e.g. through the posterior mean or mode, but also give a measure of the remaining uncertainty through the spread of the posterior.
Several papers have explored nonparametric Bayesian approaches in this setting. An early reference is Møller et al. (1998), who study log-Gaussian priors. Gugushvili and Spreij (2013) recently considered Gaussian processes combined with different, non-smooth link functions. Kernel mixtures priors are considered in Kottas and Sansó (2007). Spline-based priors are used in DiMatteo et al. (2001) and Belitser et al. (2013).
The present study is motivated by a method that is not covered by earlier theoretical papers, namely the method of Adams et al. (2009). These authors presented the first approach that is also computationally fully nonparametric in the sense that it does not involve potentially inaccurate finitedimensional approximations. The method involves a prior on the intensity that is a random multiple of a transformed Gaussian process (GP). Both the hyper parameters of the GP and the multiplicative constant are endowed with priors as well, resulting in a hierarchical Bayes procedure (details in Section 2.3). Simulation experiments and real data examples in Adams et al. (2009) show that the method can give very satisfactory results.
The aim of this paper is to advance the theoretical understanding of the method of Adams et al. (2009), which they termed "Sigmoidal Gaussian Cox Process" (SGCP). It is by now well known both from theory and practice that nonparametric Bayesian methods need to be tuned very carefully to produce good results. An unfortunate choice of the prior or incorrectly tuned hyper parameters can easily result in procedures that give misleading results or that make sub-optimal use of the information in the training data. See for instance the by now classical reference Diaconis and Freedman (1986), or the more recent paper van der Vaart and van Zanten (2011) and the references therein.
A challenge in this problem (and in nonparametric function learning in general) is to devise a procedure that avoids overfitting and underfitting. The difficulty is that the appropriate degree of "smoothing" depends on the (unknown) regularity of the intensity function that produces the data. Indeed, intuitively it is clear that if the function is very smooth then to learn the intensity at a certain location we can borrow more information from neighborrowing points than if it is very rough. Ideally we want to have a procedure that automatically uses the appropriate degree of smoothing, i.e. that adapts to regularity.
To address this issue theoretically it is common to take an asymptotic point of view. Specifically, we assume that we have n independent sets of training data, produced by Poisson processes on the d-dimensional domain S = [0, 1] d (say), with the same intensity function λ 0 : S → [0, ∞). We aim to construct the learning procedure such that we achieve an optimal learning rate, irrespective of the regularity level of the intensity. In the problem at hand it is known that if λ 0 has regularity β > 0, then the best (minimax) rate that any procedure can achieve is of the order n −β/(d+2β) , see, e.g. Kutoyants (1998) or Reynaud-Bouret (2003). Note that the smoothness degree is unknown to us, so we can not use it in the construction of the procedure, but still we want that the posterior contracts around λ 0 at the rate n −β/(d+2β) , as n → ∞, if λ 0 is β-smooth. In this paper we prove that with appropriate priors on the hyper parameters, the SGCP approach of Adams et al. (2009) attains this optimal rate (up to a logarithmic factor). It does so for every regularity level β > 0, so it is fully rate-adaptive.
Technically the paper uses the mathematical framework for studying contraction rates for Gaussian and conditionally Gaussian priors as developed in van der Vaart and van Zanten (2008a) and van der Vaart and van Zanten (2009).
We also use an extended version of a general result for Bayesian inference for 1-dimensional Poisson processes from Belitser et al. (2013). On a general level the line of reasoning is similar to that of van der Vaart and van Zanten (2009). However, due to the presence of a link function and a random multiplicative constant in the SGCP model (see Section 2 ahead) the results of the latter paper do not apply in the present setting and additional mathematical arguments are required to prove the desired results.
The paper is organised as follows. In Section 2 we describe the Poisson process observation model and the SGCP prior model, which together determine a full hierarchical Bayesian model. The main result about the performance of the SGCP approach is presented and discussed in Section 3. Mathematical proofs are given in Section 4. In Section 5 we make some concluding remarks.

The SGCP model 2.1 Observation model
We assume we observe n independent copies of an inhomogeneous Poisson process on the d-dimensional unit cube S = [0, 1] d (adaptation to other domains is straightforward). We denote these observed data by N 1 , . . . , N n .
The object of interest is the underlying intensity function. This is a (integrable) function λ : [0, 1] d → [0, ∞) with the property that given λ, every N j is a random counting measure on [0, 1] d such that N j (A) and N j (B) are independent if the sets A, B ⊂ [0, 1] d are disjoint and the number of points N j (B) falling in the set B has a Poisson distribution with mean B λ(s) ds.
If we want to stress that the probabilities and expectations involving the observations N j depend on λ, we use the notations P λ and E λ , respectively. We note that instead of considering observations from n independent Poisson processes with intensity λ, one could equivalently consider observations from a single Poisson process with intensity nλ.

Prior model
The SGCP model introduced in Adams et al. (2009) postulates a-priori that the intensity function λ is of the form where λ * > 0 is an upper bound on λ, g is a GP indexed by S and σ is the sigmoid, or logistic function on the real line, defined by σ( In the computational section of Adams et al. (2009) g is modeled as a GP with squared exponential covariance kernel and zero mean, with a prior on the length scale parameter. The hyper parameter λ * is endowed with an independent gamma prior.
In the mathematical results presented in this paper we allow a bit more flexibility in the choice of the covariance kernel of the GP, the link function σ and the priors on the hyper parameters. We assume that g is a zero-mean, homogenous GP with covariance kernel given in spectral form by where ℓ > 0 is an (inverse) length scale parameter and µ is a spectral density on R d such that the map a → µ(aξ) on (0, ∞) is decreasing for every ξ ∈ R d and that satisfies for some δ > 0 (the Euclidean inner product and norm are denoted by ·, · and · , respectively). Note that, in particular, the centered Gaussian spectral density satisfies this condition and corresponds to the squared exponential kernel Eg(s)g(t) = e −ℓ 2 t−s 2 .
We endow the length scale parameter ℓ with a prior with density p ℓ on [0, ∞), for which we assume the bounds, for positive constants C 1 , D 1 , C 2 , D 2 , nonnegative constants p, q, and every sufficiently large x > 0, This condition is, for instance, satisfied if ℓ d has a gamma distribution, which is a common choice in practice. On the upper bound λ * we put a prior satisfying an exponential tail bound. Specifically, we use a positive, continuous prior density for all λ 0 > 0. Note that this condition is fulfilled if we place a gamma prior on λ * . Finally, we use a strictly increasing, infinitely smooth link function σ : R → (0, 1) in (2.1) that satisfies for all x, y ∈ R. This condition is in particular fulfilled for the sigmoid function employed by Adams et al. (2009). It holds for other link functions as well, for instance for the cdf of the standard normal distribution.
The main results of the paper concern the posterior distribution of the intensity function λ, i.e. the conditional λ | N 1 , . . . , N n . Throughout we will denote the prior on λ by Π and the posterior by Π(· | N 1 , . . . , N n ). In this setting Bayes' formula asserts that where the likelihood is given by (see, e.g., Kutoyants (1998)).

Main result
Consider the prior and observations model described in the preceding section and let Π(· | N 1 , . . . , N n ) be the corresponding posterior distribution of the intensity function λ.
The following theorem describes how quickly the posterior distribution contracts around the true intensity λ 0 that generates the data. The rate of contraction depends on the smoothness level of λ 0 . This is quantified by assuming that λ 0 belongs to the Hölder space C β [0, 1] d for β > 0. By definition a function on [0, 1] d belongs to this space if it has partial derivatives up to the order ⌊β⌋ and if the ⌊β⌋th order partial derivatives are all Hölder continuous of the order β − ⌊β⌋. Here ⌊β⌋ denotes the greatest integer strictly smaller than β. The rate of contraction is measured in the L 2 -distance between the square root of intensities. This is the natural statistical metric in this problem, as it can be shown that in this setting the Hellinger distance between the models with intensity functions λ 1 and λ 2 is equivalent to min{ √ λ 1 − √ λ 2 2 , 1} (see Belitser et al. (2013)). Here f 2 denotes the L 2 -norm of a function on S = [0, 1] d , i.e. f 2 2 = S f 2 (s) ds. Theorem 3.1. Suppose that λ 0 ∈ C β [0, 1] d for some β > 0 and that λ 0 is strictly positive. Then for all sufficiently large M > 0, as n → ∞, for some ρ > 0.
The theorem asserts that if the intensity λ 0 that generates the data is β-smooth, then, asymptotically, all the posterior mass is concentrated in (Hellinger) balls around λ 0 with a radius that is up to a logarithmic factor of the optimal order n −β/(d+2β) . Since the procedure does not use the knowledge of the smoothness level β, this indeed shows that the method is rate-adaptive, i.e. the rate of convergence adapts automatically to the degree of smoothness of the true intensity. Let us mention once again that the conditions of the theorem are in particular fulfilled if in (2.1), λ * is taken gamma, σ is the sigmoid (logistic) function, and g is a squared exponential GP with length scale ℓ, with ℓ d a gamma variable.

Proof of Theorem 3.1
To prove the theorem we employ an extended version of a result from Belitser et al. (2013) that gives sufficient conditions for having (3.1) in the case d = 1, cf. their Theorem 1. Adaptation to the case of a general d ∈ N is straightforward. To state the result we need some (standard) notation and terminology. For a set of positive functions F we write F c for its complement and For ε > 0 and a norm · on F, let N (ε, F, || · ||) be the minimal number of balls of radius ε with respect to norm · needed to cover F. The uniform norm f ∞ of a function f on S is defined, as usual, as f ∞ = sup s∈S |f (s)|. The space of continuous function on S is denoted by C(S).
Let Π now be a general prior on the intensity function λ and let Π(· | N 1 , . . . , N n ) be the corresponding posterior (2.7).
Theorem 4.1. Assume that λ 0 is bounded away from 0. Suppose that for positive sequences δ n , δ n → ∞ such that n(δ n ∧ δ n ) 2 → ∞ 1 as n → ∞ and constants c 1 , c 2 > 0, it holds that for all L > 1, there exist subsets F n ⊂ C(S) and a constant c 3 such that Then for ε n = δ n ∨ δ n and all sufficiently large M > 0, as n → ∞.
for constants c 1 , c 2 > 0 and δ n as defined above.
The link function σ is strictly increasing and smooth, hence it has a smooth inverse σ −1 : (0, 1) → R. Define the function w 0 on S by so that λ 0 = 2 λ 0 ∞ σ(w 0 ). Since the function λ 0 positive and continuous on the compact set S it is bounded away from 0 on S, say λ 0 a > 0. It follows that λ 0 (s)/2 λ 0 ∞ varies in the compact interval [a/2||λ 0 || ∞ , 1/2] as s varies in S, whence w 0 inherits the smoothness of λ 0 , i.e. w 0 ∈ C β (S). Now observe that for ε > 0, Since λ * has a positive, continuous density the first factor on the right is bounded from below by a constant times ε. Since the function √ σ is Lipschitz by assumption, the second factor is bounded from below by P( g − w 0 ∞ cε) for a constant c > 0. By Theorem 3.1 in van der Vaart and van Zanten (2009) we have the lower bound with δ n as specified above. The proof of (4.5) is now easily completed.

Construction of sieves
Let H ℓ be the RKHS of the GP g with covariance (2.2) and let H ℓ 1 be its unit ball (see van der Vaart and van Zanten (2008b) for background on these notions). Let B 1 be the unit ball in C[0, 1] d relative to the uniform norm. Define and λ n , M n , γ n , r n and ε n are sequences to be determined later. In the next two subsections we study the metric entropy of the sieves F n and the prior mass of their complements.
(4.6) By substitutingη n = ε n √ λ n we get that for some constants K 1 and K 2 , log N (η n , F n , · 2 ) K 1 r d n log So we conclude that we have the entropy bound log N (η n , F n , · 2 ) nη 2 n for sequences λ n , M n , r n andη n satisfying (4.7) and

Remaining mass
By conditioning we have By (2.5) the second term is bounded by a constant times exp(−c 0 λ κ n ). For the first term, note that for λ λ n we have hence P(λσ(g) ∈ F n ) P(g ∈ G n ). From (5.3) in van der Vaart and van Zanten (2009) we obtain the bound P(g ∈ G n ) K 3 r p−d+1 n e −D 2 r d n log q rn log q r n + e −M 2 n /8 , for some K 3 > 0, ε n < ε 0 for a small constant ε 0 > 0, and M n , r n and ε n satisfying M 2 n > 16K 4 r d n (log(r n /ε n )) 1+d , r n > 1, (4.9) where K 4 is some large constant. It follows that P(g ∈ G n ) is bounded above by a multiple of exp (−Ln η 2 n ) for a given constant L and η n = λ n ε n , provided M n , r n , γ n and ε n satisfy (4.9) and D 2 r d n log q r n 2Ln η 2 n , r p−d+1 n e Ln η 2 n , M 2 n 8Ln η 2 n . (4.10) Note that in terms ofη n , (4.9) can be rewritten as M 2 n > 16K 4 r d n (log(r n λ n / η n )) 1+d , r n > 1. (4.11) We conclude that if (4.11),(4.10) holds and c 0 λ κ n > Ln η 2 n , (4.12) then P(λ * σ(g ∈ F n )) e −Ln η 2 n .

Completion of the proof
In the view of the preceding it only remains to show that η n ,η n , r n , M n > 1 and λ n can be chosen such that relations (4.7), (4.8), (4.10), (4.11) and (4.12) hold.
One can see that it is true for η n = δ n andη n = δ n described in the theorem, with r n , M n , λ n as follows:

Concluding remarks
We have shown that the SGCP approach to learning intensity functions proposed by Adams et al. (2009) enjoys very favorable theoretical properties, provided the priors on the hyper parameters are chosen appropriately. The result shows there is some flexibility in the construction of the prior. The squared exponential GP may be replaced by other smooth stationary processes, other link functions may be chosen, and there is also a little room in the choice of the priors on the length scale and the multiplicative parameter. This flexibility is limited, however, and although our result only gives upper bounds on the contraction rate, results like those of Castillo (2008) and van der Vaart and van Zanten (2011) lead us to believe that one might get sub-optimal performance when deviating too much from the conditions that we have imposed.
We expect that a number of generalizations of our results are be possible. For instance, it should be possible to obtain generalizations to anisotropic smoothness classes and priors as considered in Bhattacharya et al. (2014), and classes of analytic functions as studied in van der Vaart and van Zanten (2009). These generalizations take considerable additional technical work however and are therefore not worked out in this paper. We believe they would not change the general message of the paper.