© 2018 IEEE. The Non-homogeneous Poisson process is a point process with time-varying intensity across its domain, the use of which arises in numerous areas in signal processing and machine learning. However, applications are largely limited by the intractable likelihood function and the high computational cost of existing inference schemes. We present a sequential inference framework that utilises generative Poisson data and sequential Markov Chain Monte Carlo (SMCMC) algorithm to enable online inference in various applications. The proposed model is compared to competing methods on synthetic datasets and tested with real-world financial data

Godsill, SJ

Li, C

English

Apollo (Cambridge)

SEQUENTIAL INFERENCE METHODS FOR NON-HOMOGENEOUS POISSON PROCESSESWITH STATE-SPACE PRIORChenhao LiUniversity of CambridgeDepartment of EngineeringCambridge, UKSimon J. GodsillUniversity of CambridgeDepartment of EngineeringCambridge, UKABSTRACTThe Non-homogeneous Poisson process is a point processwith time-varying intensity across its domain, the use ofwhich arises in numerous areas in signal processing and ma-chine learning. However, applications are largely limited bythe intractable likelihood function and the high computationalcost of existing inference schemes. We present a sequentialinference framework that utilises generative Poisson data andsequential Markov Chain Monte Carlo (SMCMC) algorithmto enable online inference in various applications. The pro-posed model is compared to competing methods on syntheticdatasets and tested with real-world financial data.Index Terms— Bayesian inference, non-homogeneousPoisson process, state-space model, sequential-MCMC.1. INTRODUCTIONSince the original development of the Poisson process, themodel has been widely applied to various applications withpoint data in both temporal and spatial settings, such as neu-ronal spike trains [1], emissions from radioactive sourcesand traders’ operations in an open limit order book market[2]. Poisson processes have provided the possibility for de-scribing the intensity of occurrences, which is useful in boththe study of behaviour of the stochastic process itself and inthe prediction of future occurrences. The non-homogeneousPoisson process (NHPP), an important variant of the standardPoisson process, allows the intensity of the process to varywith time, giving more general and realistic modelling of thedata.Intensity inference for the NHPP has become an impor-tant research topic. The classic approach developed in earlyyears [3] uses kernel densities to construct a non-parametricestimator with empirical choice of bandwidth. The authorsof [4] presents the first tractable approach for performingBayesian inference on the Gaussian Cox process with a sig-moid transformation function. Based on these ideas, themodel was extended in [5] to include variational samplingof some hyperparameters and parallel inference for multiplecorrelated processes. Both [4] and [6] scale poorly with thenumber of data due to the high computational complexityof the Gaussian process prior. Inspired by sparse Gaussianprocess models, [5] uses inducing points to perform tractablevariational inference with the standard likelihood function,which achieves reduced complexity.In this paper, we propose a new variant of the Cox pro-cess, in which the intensity function is a transformation of thestate-space model (SSM). We then explore the SMCMC al-gorithm [7][8] as the inference method under our formulatedsequential framework. Furthermore, we extend the standardSMCMC algorithm with a mixture sampling scheme andsequential batch scheme to improve both accuracy and com-putational efficiency. Results obtained from both syntheticdatasets and a real dataset are shown and discussed in Section4 with an analysis on hyperparameter settings.2. THE MODELIn this section, we briefly describe the NHPP before review-ing the model in [4] as the Sigmoidal Gaussian Cox Pro-cess (SGCP), which achieves tractable Bayesian inference viathinning from a generative prior. Later we present our modelthat allows for the first time sequential online inference of theintensity function, by integrating the SSM and SMCMC algo-rithm into the generative framework of the SGCP model.2.1. The Non-homogeneous Poisson ProcessFor a domain S = RD of arbitrary dimension D, we candefine a non-homogeneous Poisson process with an intensityfunction λ(s), s ∈ S, for which the counting measure N(T )evaluated over a subregion T ⊂ S follows a Poisson distribu-tion with parameter λT =∫T λ(s) ds. Moreover, the numberof events in any disjoint subsets {Ti}i ⊂ S are independentrandom variables. With this definition, we can obtain the like-lihood function [9] of such a model given a set of K inputevents denoted as {sk}Kk=1 in a region T as:p({sk}Kk=1 |λ(s), T ) = exp{−∫Tds λ(s)} K∏k=1λ(sk) (1)2.2. Thinning & Tractable Joint DistributionInference for the posterior intensity is clearly intractable us-ing the likelihood in (1). However, a thinning procedure pro-vides a tractable approach for simulating a NHPP from a par-ticular intensity function λ(s) [10], and this idea forms thebasis of the SGCP inference method presented by [4]. A se-quence of events {sn}Nn=1 is generated from a homogeneousPoisson process (HPP) with intensity λ∗ ≥ λ(s), s ∈ T ,which is an upper bound on the desired NHPP intensity func-tion λ(s). Points from the desired NHPP can then be gener-ated unbiasedly by thinning the homogeneous Poisson pointsindependently with probability λ(sn)/λ∗. This constructivegeneration process leads to a tractable augmented likelihoodfunction. If we now include a mapping from a real-valuedstochastic process {g(t), t ∈ T } to the intensity function,for example the scaled sigmoidal function λ(s) = λ∗σ(g(s)),then a prior may be included over the latent values of {g(t)}.Define in ∈ {0, 1} as an indicator associated with each Pois-son event, taking value 0 for a event selected by the thinningprocess, and 1 for a ‘latent’ event rejected by the thinningprocess. A tractable joint likelihood is now obtained as:p({sn}Nn=1,g1:N , {in}Nn=1 |λ∗, T )=(λ∗)N e−λ∗|T | p(g1:N |{sn}Nn=1)N∏n=1σ{(−1)ing(sn)}(2)where {sn}Nn=1 is the time-ordered list of homogeneousPoisson points and g1:N is the vector of corresponding g(sn)stochastic process values, with p(g1:N ...) its prior probability.In [4] a Gaussian process prior [11] is adopted for p(g1:N ...)and batch-based MCMC methods are applied for inference.In our new model however a continuous-time SSM is adoptedand sequential Monte Carlo inference methods are provided,which allow for on-line inference as new data points arrive.2.3. New ModelDespite the tractability achieved in the SGCP method [4],its application is limited by its computational load, whichdemands O(N3) for the Gaussian process [11], and scalespoorly with the value of λ∗ (choice of λ∗ has also beendemonstrated to impact heavily the accuracy of the method).In order to alleviate these limitations, we propose a newmodel under the tractable framework which allows efficientsequential Bayesian inference for the intensity function. Wereplace the fully correlated Gaussian process prior with acontinuous-time SSM. This structure has the benefit thatthe prior is readily computed for arbitrary sets of Poissonpoints {sn} and also has a sequential Markovian propertythat avoids the O(N3) computational complexity and aidssequential inference formulations. We employ a SSM ofthe form dg(t) = Ag(t)dt + hdW (t) where {W (t)} isthe Wiener process. Under this framework transition den-sities are readily computed to be Gaussian and Markovian,p(g(Q)|g(P )) = N (µ(Q,P ), C(Q,P )) for Q > P , see e.g.[12], and hence the joint prior p(g1:N ) is readily computedby the probability chain rule.While any Gaussian SSM could be applied in our frame-work, we adopt a Langevin dynamics model similar to thatin [12] in which gt = [g1,t g2,t]T contains a stochastic trendterm g2,t that is integrated to give the value g1,t which is in-put to the sigmoidal function of the SGCP framework. Inthis model A =[0 10 θ]and h = [0 σ]T , where θ ≥ 0 isa mean-reverting term on the trend and σ is a scale param-eter for the noise W (t). The joint prior under this model isthen computed as in [12] and denoted p(g1:N |{sn}Nn=1) =LD(g1:N |{sn}Nn=1).3. SEQUENTIAL INFERENCEInference can be made sequential by inputting short batchesof data, delineated by times tk, i = 0, 1, .... The kth timeinterval is Ek = (tk−1, tk] ⊂ T , and for example we can useregularly spaced intervals tk = kδT , or we can space accord-ing to the arrival times of observed points {sn; in = 0}. Wefurther define the notation {s,g, i}Ek = {sn,gn, in; sn ∈Ek} as the locations, the state vectors and the indicators cor-responding to all events in the interval Ek. We can hence writethe joint distribution recursion as:p({s,g, i}E1:k |λ∗, E1:k) = p({s,g, i}E1:k−1 |λ∗, E1:k−1)× p({s,g, i}Ek |λ∗, {s,g, i}Ek−1 , Ek)(3)where the conditional propagation can be factorised based on(2) as:p({s,g, i}Ek |λ∗, {s,g, i}Ek−1 , Ek) = (λ∗)Nke−λ∗|Ek|× LD({g}Ek |{g}Ek−1)×∏n:sn∈Ekσ{(−1)ing1,n} (4)where Nk = |{n; sn ∈ Ek}| is the total number of eventsin Ek. A standard (variable rate) particle filter can be usedfor inference but is not efficient as the proposal in (4) in-volves multiple latent variables in a single propagation. Weaddress this high-dimensional proposal with the SMCMCalgorithm where samples are obtained with local or globalMCMC moves followed by a Metropolis-Hastings (MH)accept-reject step as bias correction [13][14]. Furthermore,a mixture sampling procedure is adopted: at each MCMCiteration, a decision is made on performing either a joint MHproposal step with probability PJ or a sequence of individualrefinement Metropolis-within-Gibbs transitions with proba-bility 1−PJ . This scheme provides trade-offs between theefficiency and the accuracy of the inference.3.1. Joint ProposalA joint MH kernel consists of a discrete uniform draw of{g, s}(p)Ek−1 from the particle collection of Ek−1 followed bythe proposals of {g, s}Ek conditioned on the drawn particlep. More specifically, the latter is achieved by three samplersin combination: 1) the total number of thinned events M̃ in Ekfrom a Poisson distribution; 2) the locations of thinned events{s̃m}M̃m=1 uniformly from Ek; and 3) state vectors {g}Ek ofall events in Ek from the LD prior conditioned on all eventlocations {s}Ek and the drawn particle {g, s}(p)Ek−1 . This jointproposal has density:qJ =Poi(M̃ |λ∗, Ek)|Ek|M̃LD({g}Ek |{g}Ek−1) (5)Incorporating (3) and (4), we obtain the MH acceptance ratiofor the jth MCMC iteration as:ρJ =qj−1J (λ∗)N∗kLD({g}∗Ek |{g}∗Ek−1)∏n σ{(−1)i∗ng∗1,n}q∗J(λ∗)Nj−1k LD({g}j−1Ek |{g}j−1Ek−1)∏n σ{(−1)ij−1n gj−11,n}(6)where the superscript ∗ indicates the variables proposed incurrent iteration and j−1 indicates the accepted variables fromlast iteration. Nburn iterations are run before before include theaccepted particle into new particle set to ensure convergence.Tuning proposals with domain knowledge can certainly im-prove convergence rate and performance.3.2. Refinement Metropolis-within-GibbsA Joint proposal can sometimes result in low acceptance ra-tios and ultimately low particle diversity. We hence use threelocal conditional samplers to individually move the threetypes of latent variables. Taking advantages from the SGCPmodel, we use a similar construction as presented in [4]: aMH move to perturb the number of thinned events M̃ in Ek;a MH move to change the locations of the thinned events{s̃m}M̃m=1 conditioned on M̃ ; and a Metropolis-Adjusted-Langevin-Algorithm (MALA) for state vectors {g}Ek tomake efficient use of the gradient information available. Therefinement moves improve inference accuracy at a cost ofextra computation.3.3. Sequential Batch SchemeWe defined earlier the interval Ek. If the times of the batchintervals are fixed at the observed event times then this canbe non-ideal for inference as it decouples the local locationinformation among non-neighbouring points. Hence, in thisTable 1. Numerical results for models. Bold is the best.S-LD KDE SGCPλ1(s)mse 0.0257 0.129 0.0704L(p) 1.825 – -9.440Time (s) 15.86 0.01 60.23λ2(s)mse 0.6531 0.8599 1.5257L(p) -248.1 – -326.6Time (s) 60.05 0.05 1326.28(1) Data and model fits for λ1(s)(2) Data and model fits for λ2(s)Fig. 1. Model performances on synthetic datasetssection we experiment with regular sized batches of durationδT . With an appropriate choice of the batch size, the schemerecovers the temporal correlation among points in the samebatch and thus improves the sequential inference accuracy.Furthermore by placing a Gamma prior over the maximumintensity λ∗, we obtain a posterior Gibbs update of the lo-cal λ∗k for each batch, which provides considerable compu-tational saving compared to a global λ∗: αpost = α + Nk,βpost = β + |Ek|.4. RESULTS AND DISCUSSIONSIn this section, we present three empirical analyses of oursequential-Langevin (S-LD) model. We use synthetic datasetswith ground truth λ(s) to compare the relative performance ofS-LD, SGCP and a baseline kernel density estimation (KDE)method [3]. We then apply our model to a financial datasetwith high frequency input of events. Finally, we examine theeffects of hyperparameters on the S-LD model performance.Fig. 2. Results of the S-LD model on limit order book data.The top plot shows the inferred limit order arrival intensity atprice $1.12960. The bottom plot shows the market mid-priceand the selected price level for the same duration4.1. Synthetic DataTwo sets of one-dimensional data are created with the follow-ing intensity functions:1. A sum of an exponential and a Gaussian bump:λ1(s) = 2 exp{−s/15} + exp{−((s − 25)/10)2} onthe interval [0, 50] with 55 events.2. A doubly-stochastic process with λ2(s) governed byLangevin dynamics with parameters θ=−0.5, σ=0.5on interval [0, 100] with 156 events.We compared the S-LD model to both the SGCP model witha square-exponential covariance function and the KDE ap-proach with Gaussian smoothing kernel. Fig.1 shows thegraphical results of the three models and Table.1 quantita-tively reports the performance averaged across 10 trials interms of the mean squared error (MSE), the (log) likelihoodof truth under Gaussian assumption and computational time.Visually and numerically, the S-LD outperforms the othertwo methods in both MSE and likelihood. Computational-wise, our model gives reasonable computational costs evenwith high number of input points whilst the SGCP modelis significantly more expensive with λ2(s). The KDE givesthe best computational speed but coarse estimations of theintensity overall.4.2. Application to Order Book DataWe ran the S-LD model on a small set of limit order book(LOB) data taken from the EUR-USD FOREX market on the2nd of September 20151. The data describes the limit orderarrivals at a fixed price level of $1.12960 for a duration of5 minutes (19:35–19:40) from one of the busiest hour of theday. The specific interval and price is picked intentionallyto show several transitions between bid and ask sides. Fig.2shows the outcome from the S-LD model. Note the SGCPmodel is far too expensive on this dataset.1The authors would like to thank Cambridge Capital Management forproviding the datasets for these experiments.0.05 0.2 0.4 0.6 0.8 1.0Joint Proposal Prob. PJ0.00.10.20.30.40.50.6MSEMSE0.010.020.030.040.050.060.070.0Time (s)MSE and Time for different settings of PJComp. Cost(1) Varying PJ2.0 4.0 6.0 8.0 10.0No. of Batches0.00.050.10.150.20.250.30.350.40.45MSEMSE30.035.040.045.050.055.060.065.070.075.0Time (s)MSE and Time for different number of batchesComp. Cost(2) Varying the number of batchesFig. 3. The S-LD model performances under different settingsof hyperparameters under synthetic dataset λ1(s)4.3. Hyperparameter SettingsThis analysis focuses on the scheme-related hyperparametersinstead of the hyperparameters of the SSM as the latter canbe made adaptive with extension of a variational structure[15]. The two plots in Fig.3 show how computational timeand MSE change with different settings of PJ and the num-ber of batches respectively while other hyperparameters areheld fixed. The inference accuracy is improved considerablyupon the use of the refinement samplers. Reducing PJ lin-early increases the computation time and eventually raises theMSE due to the lack of sampling from previous particle col-lection. The sequential batch scheme improves both accuracyand efficiency of the inference. However, too many batchesgive similar outcomes as the pointwise propagation scheme.The sudden jump between PJ=1.0 and 0.9 indicates that therefinement procedure significantly improves the acceptanceratio by proposing samples in a Gibbs manner.5. SUMMARYWe have introduced a novel sequential method for inferenceabout the intensity function in a NHPP. By avoiding the in-tractability with generative prior and latent variables, ourmodel has demonstrated improved accuracy and efficiencycompared to the KDE and SGCP approaches. The sequentialframework utilising SMCMC algorithm allows our modelto stand in contrast to other approaches for the NHPP inthat it provides online inference and the possibility of auto-adaptation of hyperparameters to input data, which can beachieved with a variational structure. The sequential batchscheme further improves the model performance by restor-ing the local location information, which is crucial in theinference of NHPP but usually ignored in standard sequentialinference methods.6. REFERENCES[1] D. Perkel, G. Gerstein, and G. Moore, “Neuronal spiketrains and stochastic point processes: II. Simultaneousspike trains,” Biophysical journal, vol. 7, no. 4, pp. 419–440, 1967.[2] M. Avellaneda and S. Stoikov, “High-frequency tradingin a limit order book,” Quantitative Finance, vol. 8, no.3, pp. 217–224, 2008.[3] P. Diggle, “A kernel method for smoothing point processdata,” Applied statistics, pp. 138–147, 1985.[4] R. Adams, I. Murray, and D. MacKay, “Tractable non-parametric Bayesian inference in Poisson processes withGaussian process intensities,” in Proceedings of the 26thAnnual International Conference on Machine Learning.ACM, 2009, pp. 9–16.[5] C. Lloyd, T. Gunter, M. Osborne, and S. Roberts, “Vari-ational inference for Gaussian process modulated Pois-son processes,” in International Conference on MachineLearning, 2015, pp. 1814–1822.[6] T. Gunter, C. Lloyd, M. Osborne, and S. Roberts, “Ef-ficient Bayesian nonparametric modelling of structuredpoint processes,” arXiv preprint arXiv:1407.6949, 2014.[7] S. K. Pang, J. Li, and S. Godsill, “Detection andtracking of coordinated groups,” IEEE Transactions onAerospace and Electronic Systems, vol. 47, no. 1, pp.472–502, 2011.[8] A. Golightly and D. Wilkinson, “Bayesian sequential in-ference for nonlinear multivariate diffusions,” Statisticsand Computing, vol. 16, no. 4, pp. 323–338, 2006.[9] D. Cox and V. Isham, Point processes, vol. 12, CRCPress, 1980.[10] P. Lewis and G. Shedler, “Simulation of nonhomoge-neous Poisson processes by thinning,” Naval ResearchLogistics (NRL), vol. 26, no. 3, pp. 403–413, 1979.[11] C. Rasmussen and C. Williams, Gaussian processes formachine learning, vol. 1, MIT press Cambridge, 2006.[12] H. Christensen, J. Murphy, and S. Godsill, “Forecastinghigh-frequency futures returns using online Langevindynamics,” IEEE Journal of Selected Topics in SignalProcessing, vol. 6, no. 4, pp. 366–380, 2012.[13] S. K. Pang, S. Godsill, J. Li, F. Septier, and S. Hill, Se-quential inference for dynamically evolving groups ofobjects, chapter 12, pp. 245–276, Bayesian Time Se-ries Models. 2011.[14] F. Septier and G. Peters, “Langevin and Hamiltonianbased sequential MCMC for efficient Bayesian filteringin high-dimensional spaces,” IEEE Journal of SelectedTopics in Signal Processing, vol. 10, no. 2, pp. 312–327,2016.[15] Z. Ghahramani and G. Hinton, “Variational learningfor switching state-space models,” Neural computation,vol. 12, no. 4, pp. 831–864, 2000.

Sequential inference methods for non-homogeneous poisson processes with state-space prior

Crossref

Sequential Inference Methods for Non-Homogeneous Poisson Processes with State-Space Prior

https://www.repository.cam.ac.uk/bitstream/handle/1810/317063/ICASSP_published.pdf?sequence=1&isAllowed=y

Sequential inference methods for non-homogeneous poisson processes with state-space prior

Abstract

Similar works

Full text

Available Versions

Apollo (Cambridge)

Crossref