Classification and clustering for observations of event time data using non-homogeneous Poisson process models

TL;DR

This paper develops methods for classifying and clustering event time data modeled by non-homogeneous Poisson processes, using spline-based rate function estimation and maximum likelihood techniques, validated on synthetic and real data.

Contribution

It introduces a spline-based modeling approach for NHPP rate functions and applies maximum likelihood and EM algorithms for classification and clustering tasks.

Findings

Methods perform well on synthetic data.

Effective classification and clustering on real-world data.

Code available for reproducibility.

Abstract

Data of the form of event times arise in various applications. A simple model for such data is a non-homogeneous Poisson process (NHPP) which is specified by a rate function that depends on time. We consider the problem of having access to multiple independent observations of event time data, observed on a common interval, from which we wish to classify or cluster the observations according to their rate functions. Each rate function is unknown but assumed to belong to a finite number of rate functions each defining a distinct class. We model the rate functions using a spline basis expansion, the coefficients of which need to be estimated from data. The classification approach consists of using training data for which the class membership is known, to calculate maximum likelihood estimates of the coefficients for each group, then assigning test observations to a group by a maximum likelihood criterion. For clustering, by analogy to the Gaussian mixture model approach for Euclidean data, we consider mixtures of NHPP and use the expectation-maximisation algorithm to estimate the coefficients of the rate functions for the component models and group membership probabilities for each observation. The classification and clustering approaches perform well on both synthetic and real-world data sets. Code associated with this paper is available at https://github.com/duncan-barrack/NHPP .

Equations

$$l(\lambda (t)|\mathbf{\chi })=-M(t)+\sum _{i=1}^{m_{\chi }}\text{mbox}ln(\lambda ({\chi }_i)),$$

$$\lambda (t)=\sum _{m=1}^{n_b}{c}_m{B}_m(t).$$

$$f(\mathbf{c})=-l(\lambda (t)|\mathbf{\chi })=M(t)-\sum _{i=1}^{m_{\chi }}\text{mbox}ln\left(\sum _{m=1}^{n_b}{c}_m{B}_m({\chi }_i)\right),$$

$$g(\mathbf{c})=-\lambda (t)=-\sum _{m=1}^{n_b}{c}_m{B}_m(t)\le 0.$$

$$▽f$$

$$\mathbf{H}=\left[\begin{array}{ccc}{H}_{1,1}(f)& \dots & H_{1,n_b}(f)\\ ⋮& \ddots & ⋮\\ H_{n_b,1}(f)& \dots & H_{n_b,n_b}(f)\end{array}\right],$$

$$\left[\begin{array}{ccc}\frac{B_1({\chi }_1)}{{\sum }_{m=1}^{n_b}{c}_m{B}_m({\chi }_1)}& \dots & \frac{B_{n_b}({\chi }_1)}{{\sum }_{m=1}^{n_b}{c}_m{B}_m({\chi }_1)}\\ ⋮& \ddots & ⋮\\ \frac{B_1({\chi }_{m_{\chi }})}{{\sum }_{m=1}^{n_b}{c}_m{B}_m({\chi }_{m_{\chi }})}& \dots & \frac{B_{n_b}({\chi }_{m_{\chi }})}{{\sum }_{m=1}^{n_b}{c}_m{B}_m({\chi }_{m_{\chi }})}\end{array}\right],$$

$$l({\lambda }_{\nu }(t)|{\mathbf{T}}_{\text{train}})=$$

$${\lambda }_{\nu }(t)=\sum _{m=1}^{n_b}{c}_{\nu ,m}{B}_m(t).$$

$$P(g_j=\nu |{\mathbf{\chi }}_j)=\frac{e^{-{\hat{M}}_{\nu }(t)}{\prod }_{i=1}^{m_j}{\hat{\lambda }}_{\nu }({\chi }_{j,i})}{{\sum }_{q=1}^k{e}^{-{\hat{M}}_q(t)}{\prod }_{i=1}^{m_j}{\hat{\lambda }}_q({\chi }_{j,i})}.$$

$$\begin{array}{rl}l(\mathbf{\theta }|\mathbf{D},\mathbf{z})& =\sum _{l=1}^n\text{mbox}ln\left(p({\mathbf{\chi }}_l,\mathbf{z}|\mathbf{\theta })\right)\\ & =\sum _{l=1}^n\text{mbox}ln\left(\sum _{q=1}^k\text{mathds}1(z_l=q){\tau }_q\mathcal{L}({\lambda }_q(t)|{\mathbf{\chi }}_l)\right),\end{array}$$

$$\begin{array}{rl}l(\mathbf{\theta }|\mathbf{D})& =\sum _{l=1}^n\text{mbox}ln\left(p({\mathbf{\chi }}_l|\mathbf{\theta })\right)\\ & =\sum _{l=1}^n\text{mbox}ln\left(\sum _{q=1}^k{\tau }_q\mathcal{L}({\lambda }_q(t)|{\mathbf{\chi }}_l)\right).\end{array}$$

$$Q(\mathbf{\theta },{\hat{\mathbf{\theta }}}^c)$$

$$=\sum _{l=1}^n\text{mbox}{E}_{\mathbf{z}|\mathbf{D},{\hat{\mathbf{\theta }}}^c}[l(\mathbf{\theta }|{\mathbf{\chi }}_l,z_l)],$$

$$=\sum _{l=1}^n\sum _{q=1}^{k}p(z_l=q|{\mathbf{x}}_{\mathbf{l}},{\hat{\mathbf{\theta }}}^c)l(\mathbf{\theta }|{\mathbf{x}}_{\mathbf{l}},z_l),$$

$$=\sum _{l=1}^n\sum _{q=1}^k{\hat{R}}_{q,l}^c\left(\text{mbox}ln({\tau }_q)+l({\lambda }_q(t)|{\mathbf{\chi }}_l))\right).$$

$$Q(\mathbf{\theta },{\hat{\mathbf{\theta }}}^c)=\sum _{l=1}^n\sum _{q=1}^k{\hat{R}}_{q,l}^c\left(\text{mbox}ln({\tau }_q)-M_q(t)+\sum _{i=1}^{m_l}\text{mbox}ln({\lambda }_q({\chi }_{l,i}))\right),$$

$$\begin{array}{rl}{\hat{R}}_{\nu ,l}^c& =p(z_l=\nu |{\mathbf{\chi }}_l,{\hat{\mathbf{\theta }}}^c)\\ & =\frac{p(z_l=\nu )p({\mathbf{\chi }}_l|z_l=\nu ,{\hat{\mathbf{\theta }}}^c)}{{\sum }_{q=1}^{k}p(z_l=q)p({\mathbf{\chi }}_l|z_l=q,{\hat{\mathbf{\theta }}}^c)}\\ & =\frac{{\hat{\tau }}_{\nu }^c{e}^{-{\hat{M}}_{\nu }^c(t)}{\prod }_{i=1}^{m_l}{\hat{\lambda }}_{\nu }^c({\chi }_{l,i})}{{\sum }_{q=1}^k{\hat{\tau }}_q^c{e}^{-{\hat{M}}_q^c(t)}{\prod }_{i=1}^{m_l}{\hat{\lambda }}_q^c({\chi }_{l,i})}.\end{array}$$

$$\mathcal{L}(\mathbf{\theta },{\hat{\mathbf{\theta }}}^c,\delta )=\sum _{l=1}^n\sum _{q=1}^k{\hat{R}}_{q,l}^c\left(\text{mbox}ln({\tau }_q)-M_q(t)+\sum _{i=1}^{m_l}\text{mbox}ln({\lambda }_q({\chi }_{l,i}))\right)+\delta \left(\sum _{q=1}^k{\tau }_q-1\right),$$

$$\frac{\partial \mathcal{L}(\mathbf{\theta },{\mathbf{\theta }}^c,\delta )}{\partial {\tau }_{\nu }}=$$

$$+\frac{\partial }{\partial {\tau }_{\nu }}\delta \left(\sum _{q=1}^k{\tau }_q-1\right).$$

$$\frac{\partial \mathcal{L}(\mathbf{\theta },{\mathbf{\theta }}^c,\delta )}{\partial {\tau }_{\nu }}$$

$$=\sum _{l=1}^n\frac{{\hat{R}}_{\nu ,l}^c}{{\tau }_{\nu }}+\delta .$$

$${\tau }_{\nu }=\sum _{l=1}^n\frac{{\hat{R}}_{\nu ,l}^c}{-\delta }.$$

$$\sum _{m=1}^{n_b}{c}_{q,m}{B}_m(t),\dots ,\sum _{m=1}^{n_b}{c}_{k,m}{B}_m(t)\ge 0,$$

$${\lambda }_1(t)$$

$$\lambda_{3}(t)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}20&\text{if $t<\frac{t_{2}}{4}$}\\ 40&\text{if $\frac{t_{2}}{4}\leq t<\frac{t_{2}}{2}$}\\ 60&\text{if $\frac{t_{2}}{3}\leq t<\frac{3t_{2}}{4}$}\\ 80&\text{if $t\geq\frac{3t_{2}}{4}$}\\ \end{array}\right.$$

$$\lambda_{4}(t)=\left\{\begin{array}[]{@{}c@{\quad}l@{}}80&\text{if $t<\frac{t_{2}}{4}$}\\ 60&\text{if $\frac{t_{2}}{4}\leq t<\frac{t_{2}}{2}$}\\ 40&\text{if $\frac{t_{2}}{3}\leq t<\frac{3t_{2}}{4}$}\\ 20&\text{if $t\geq\frac{3t_{2}}{4}$}\\ \end{array}\right.$$

Full text

1 Introduction

Much real world data is amenable to modelling as a non-homogeneous Poisson process (NHPP) (e.g. stock purchase transaction times (Engle, 2000), email traffic (Malmgren et al., 2008) as well as computer network traffic (Paxson and Floyd, 1995)). In various applications a natural unit of observation is a sample of event times observed on some common interval $\{t:0<t\leq T\}$ and a challenge is to label each observation as belonging to one of a small number of distinct classes. A particular example we consider later is the categorisation of stores by till transaction times. This can provide retailers valuable insights in to the different demands placed on different store types and help to inform various decisions from stock management to appropriate store staffing levels. Another example we consider is the categorisation of transportation hubs according to the times at which they are used. Such information helps inform the type and number of services (transportation, catering, retail and others) that should be made available at each type of hub. When a ‘training set’, containing observations with known class memberships is available, such a task is termed ‘classification’, when no labels are available for the data it is termed a ‘clustering’ problem (Hastie et al., 2009). Previous approaches to classification and clustering of event time data have focused on deriving summary statistics from the data such as the mean and variance of the inter-event times and using these as inputs in statistical learning models (Moore and Zuev, 2005, Roughan et al., 2004, McGregor et al., 2004). In addition to studies of temporal data a number of authors have considered spatial point processes. Illian et al used functional principal component analysis on second order summary statistics of spatial point process data to classify plant species (Illian et al., 2006, 2004). Closer in approach to the present work Linnett et al. (1995) used a maximum likelihood (ML) criterion to classify the textures of images modelled using homogeneous Poisson processes (HPPs).

HPPs, like those used in the work of Linnett et al. (1995) assume rates are constant for each class which is too restrictive for many applications. For example stock transaction times and human communication behaviour both exhibit periodic patterns (Engle, 2000, Malmgren et al., 2008). Computer network traffic is better modelled by non-homogeneous NHPP models rather than HPPs (Paxson and Floyd, 1995). Hence, in this work, we consider instead NHPP models, which are defined in terms of rate functions that vary with time. We model these rate functions using a spline basis expansion which imposes an assumption of smoothness on the functions and means estimating them amounts to estimating a finite number of basis coefficients. We go beyond previous work which has focused purely on rate function estimation (Alizadeh et al., 2008, Kuhl and Bhairgond, 2000, Scott et al., 1980, Zhao and Xie, 1996, Massey et al., 1996, Lee et al., 1991) by considering the tasks of classification and clustering. For classification we show that once intensity estimates for the NHPP rate functions corresponding to each class are obtained by solving a convex optmisation problem, it is a simple task to obtain the posterior class membership probabilities for test data. For clustering, where class labels are not known a priori, we assume all data to have been generated from a mixture of $k$ NHPP models. Again the rate functions for each model are represented by a combination of basis functions but his time the expectation-maximisation (EM) algorithm is used to obtain the coefficients of the basis functions as well as $k$ membership probabilities for each distinct observation of event times. Our approaches are validated on synthetic and three real data sets.

The paper is structured as follows. In Section 2 we describe our procedure for estimating NHPP rate functions. Then in Sections 3 and 4 we outline our classification and clustering methods. The results of applying these to synthetic and real data are shown in Sections 5. We conclude with a discussion in Section 6.

2 Rate function estimation

Suppose that we assume that a sample of event times $\boldsymbol{\chi}=\{\chi_{1},\chi_{2},\ldots,\chi_{m_{\chi}}\}$, observed on the interval $\{t:0<t\leq T\}$, arises from a NHPP with rate function $\lambda(t)$, then the log likelihood for the rate function is

[TABLE]

where $M(t)=\int_{0}^{t}\lambda(\tilde{t})d\tilde{t}$ (see for example Thompson (2012)).

In estimating $\lambda(t)$ in equation (1) we impose that it belongs to a class of smooth functions by assuming it is a linear combination of smooth basis functions

[TABLE]

In this paper, for the basis functions $B_{1}(t),\ldots,B_{n_{b}}(t)$, we use B-splines (De Boor, 1978) which offer control over smoothness and have finite support. We estimate the basis function coefficients $\boldsymbol{c}=\{c_{1},c_{2},\ldots,c_{n_{b}}\}$ by minimising $-l(\lambda(t)|\boldsymbol{\chi})$ i.e. by minimising

[TABLE]

subject to the constraint that the rate function $\lambda(t)$ is non-negative i.e.

[TABLE]

Minimising (3) subject to (4) is a convex optimisation problem which we show as follows. First, we derive the gradient of $f(\boldsymbol{c})$

[TABLE]

where $\bigtriangledown_{j}(f)=\frac{\partial}{\partial c_{j}}M(t)-\sum_{i=1}^{m_{\chi}}\frac{B_{i}(\chi_{i})}{\sum_{m=1}^{n_{b}}c_{m}B_{m}(\chi_{i})}$. The $n_{b}\times n_{b}$ Hessian is

[TABLE]

where $H_{s,u}(f)=\sum_{i=1}^{m_{\chi}}\frac{B_{s}(\chi_{i})B_{u}(\chi_{i})}{\left(\sum_{m=1}^{n_{b}}c_{m}B_{m}(\chi_{i})\right)^{2}}$

as $\frac{\partial^{2}}{\partial c_{s}\partial c_{u}}M(t)=\int_{t_{1}}^{t}\left(\frac{\partial^{2}}{\partial c_{s}\partial c_{u}}\sum_{m=1}^{n_{b}}c_{m}B_{m}(\tilde{t})\right)d\tilde{t}=0$ for every $s$ and $u$.

As we may express $\mathbf{H}=\mathbf{X}^{\prime}\mathbf{X}$, where the $\chi_{m_{\chi}}\times n_{b}$ matrix $\mathbf{X}=$

[TABLE]

thus $\mathbf{H}$ is positive semi-definitive as $\boldsymbol{y}^{\prime}\mathbf{H}\boldsymbol{y}=\boldsymbol{y}^{\prime}\mathbf{X}^{\prime}\mathbf{X}\boldsymbol{y}=(\boldsymbol{y}\mathbf{X})^{\prime}(\mathbf{X}\boldsymbol{y})=||\mathbf{X}\boldsymbol{y}||^{2}\geq 0$ for every non-zero column vector $\boldsymbol{y}=\{y_{1},y_{2},\ldots,y_{n_{b}}\}$ of real numbers. Therefore $f(\boldsymbol{c})$ (equation (3)) is convex in $\boldsymbol{c}$.

Constraint (4) is linear in $\boldsymbol{c}$ and hence convex in $\boldsymbol{c}$. Hence estimating the basis function coefficients of the NHPP rate function $\lambda(t)$ is a convex optimisation problem.

Furthermore, as $f(\boldsymbol{c})$ and $g(\boldsymbol{c})$ are twice differentiable, the minimisation to perform is straightforward using anyone of several numerical techniques designed for solving non-linear convex optimisation problems where the objective function and inequality constraint are both twice differentiable. For the results in this paper we used the interior point method (Byrd et al., 1999, 2000, Waltz et al., 2006).

3 Classification

Consider a set of $n_{\text{train}}$ training observations $\boldsymbol{T}_{\text{train}}=$

$\{\{\boldsymbol{\chi}_{1},g_{1}\},\{\boldsymbol{\chi}_{2},g_{2}\},\ldots,\{\boldsymbol{\chi}_{n_{\text{train}}},g_{n_{\text{train}}}\}\}$ in which $\boldsymbol{\chi}_{l}=\{\chi_{l,1},\chi_{l,2},\ldots,\chi_{l,m_{l}}\}$ is a set of of $m_{l}$ event times for the $l^{\text{th}}$ observation and $g_{l}\in{1,2,\ldots,k}$ is its corresponding class label. The classification task consists of using this training data to predict the unobserved class label for a new observation $j$ with event times $\boldsymbol{\chi}_{j}=\{\chi_{j,1},\chi_{j,2},\ldots,\chi_{j,m_{j}}\}$. Assuming each training observation is independent, the log likelihood estimate for the rate function $\lambda_{\nu}(t)$ for class $\nu$ is

[TABLE]

where the indicator function $\mathds{1}(G_{l}=\nu)=1$ if $G_{l}=\nu$ and 0 otherwise. The rate function $\lambda_{\nu}(t)$ is given by

[TABLE]

Estimates for the basis function coefficients of this expression are found via the maximum likelihood estimation outlined in Section 2.

**Classification of test data

**The posterior probability that a test observation $\boldsymbol{\chi}_{j}$ is a member of class $\nu$ can be derived from Bayes’ theorem

[TABLE]

Observations are then assigned to the class for which this probability is maximal, i.e. $\widehat{g}=\underset{\nu}{\mathrm{argmax}}\{P(g_{j}=\nu|\boldsymbol{\chi}_{j})\}$.

4 Clustering

For the unsupervised learning task of cluster analysis of observations of event time samples we assume the data to have been generated from a mixture of NHPP models. Furthermore, although class memberships are unknown, each observation is associated with a hidden latent variable which specifies which group, or class, it belongs to. Specifically, we assume the event times samples from $n$ observations $\boldsymbol{D}=\{\boldsymbol{\chi}_{1},\boldsymbol{\chi}_{2},\ldots,\boldsymbol{\chi}_{n}\}$ are generated from a mixture of $k$ NHPP models and let $\boldsymbol{z}=\{z_{1},z_{2},\ldots,z_{n}\}$ be a set of latent random variables which determine the component from which each observation originates. For the mixture the complete data log likelihood function is

[TABLE]

and corresponding incomplete data log likelihood is

[TABLE]

Here $\tau_{\nu}$ is the ‘mixing weight’ of component $\nu$, where $0\leq\tau_{\nu}\leq 1$, $\sum^{k}_{q=1}\tau_{q}=1$, and $\mathcal{L}(\lambda_{\nu}(t)|\boldsymbol{\chi}_{l})$ is the likelihood for the rate function associated with NHPP model $\nu$ given the observation $\boldsymbol{\chi}_{l}$. As before, we represent the rate functions associated with each component $\nu$ as a linear combination of basis functions, i.e. $\lambda_{\nu}(t)=\sum_{m=1}^{n_{b}}c_{\nu,m}B_{m}(t)$.

We use the EM algorithm (see Dempster et al. (1977) for details) to find estimates for the parameters $\boldsymbol{\theta}=\{\boldsymbol{\tau},\boldsymbol{C}\}$ (where $\boldsymbol{\tau}=\{\tau_{1},\tau_{2},\ldots,\tau_{k}\}$, $\boldsymbol{C}=\{\boldsymbol{c}_{1},\boldsymbol{c}_{2},\ldots,\boldsymbol{c}_{k}\}$ and $\boldsymbol{c}_{\nu}=\{c_{\nu,1},\ldots,c_{\nu,n_{b}}\}$ is the set of basis function coefficients for NHPP rate function for group $\nu$) which maximise the log likelihood of the mixture. This procedure is described below.

**Expectation step

**Given a current estimate for the model parameters $\widehat{\boldsymbol{\theta}}^{c}$ we find the so called ‘auxiliary function’ which is the expected value of the data log likelihood (LABEL:eq:lik_EM_comlete) of the mixture model with respect to the conditional distribution of $\boldsymbol{z}$ given the data $\boldsymbol{D}$.

[TABLE]

Here $\widehat{R}_{\nu,l}^{c}=p(z_{l}=\nu|\boldsymbol{\chi_{l}},\widehat{\boldsymbol{\theta}}^{c})$ is the ‘membership probability’ of observation $l$ to group $\nu$. $l(\lambda_{\nu}(t)|\boldsymbol{\chi}_{l})=-M_{\nu}(t)+\sum_{i=1}^{m_{l}}\lambda_{\nu}(\chi_{l,i})$ is the log likelihood for rate function $\nu$ for observation $l$. Plugging in this latter expression for the log likelihood into (13) we get

[TABLE]

To obtain the membership probabilities for each observation $l$ of model $\nu$ we regard the mixing weights $\boldsymbol{\tau}$ as prior probabilities of each mixture component and use Baye’s theorem to obtain

[TABLE]

**Maximisation step

**After obtaining values for the membership probabilities, we obtain subsequent estimates for the parameters $\widehat{\boldsymbol{\theta}}^{c+1}$ by finding the $\boldsymbol{\theta}$ which maximise the auxiliary function (14) subject to the constraint that the mixing weights must sum to 1, i.e. $\sum^{k}_{q=1}\tau_{q}=1$. To address this constraint we consider the Lagrange function

[TABLE]

where $\delta$ is the Lagrange multiplier. Updates for the mixing weights $\boldsymbol{\tau}$ are found by taking the partial derivative of (16) with respect to $\tau_{\nu}$, setting it to zero and solving for $\boldsymbol{\tau}$. The partial derivative of $\mathcal{L}(\boldsymbol{\theta},\boldsymbol{\theta}^{c},\delta)$ is

[TABLE]

Evaluating (17) and removing all constant terms we get

[TABLE]

Setting equation (18) to 0 and solving for $\tau_{\nu}$ we get

[TABLE]

Using the constraint that $\sum^{k}_{q=1}\tau_{q}=1$, then $-\delta=\sum_{l=1}^{n}\sum^{k}_{q=1}\widehat{R}_{q,l}^{c}=\sum_{l=1}^{n}1=n$. Hence the updated estimate for the mixing weight is $\widehat{\tau}_{\nu}^{c+1}=\sum_{l=1}^{n}\frac{\widehat{R}_{\nu,l}^{c}}{n}$.

Next we seek the $\boldsymbol{C}$ which maximise $\mathcal{L}(\boldsymbol{\theta},\boldsymbol{\theta}^{c},\delta)$ subject to the constraints that all rate functions are non-negative, i.e.

[TABLE]

for all $q=1,\ldots,k$.

This is a convex minimisation problem of the same variety as encountered in Section 2. To see this, consider the term $-f(\boldsymbol{c}_{\nu,l})$=

$-M_{q}(t)+\sum_{i=1}^{m_{l}}\mbox{ln}\left(\sum_{m=1}^{n_{b}}c_{q,m}B_{m}(\chi_{l,i})\right)$ from equation (16). $f(\boldsymbol{c}_{\nu},l)$ is identical to equation (3) which we showed is convex in the set of basis function coefficients in Section 2. As $\widehat{R}_{\nu,l}\geq 0$ for every $\nu$ and $l$, $\widehat{R}_{\nu,l}f(\boldsymbol{c}_{\nu,l})=$

$\widehat{R}_{\nu,l}\left(M_{q}(t)-\sum_{i=1}^{m_{l}}\mbox{ln}(\sum_{m=1}^{n_{b}}c_{q,m}B_{m}(\chi_{l,i}))\right)$ is also convex. Furthermore, disregarding terms which do not involve $\boldsymbol{C},-\mathcal{L}(\boldsymbol{\theta},\widehat{\boldsymbol{\theta}}^{c},\delta)=\sum_{l=1}^{n}\sum^{k}_{q=1}\widehat{R}_{q,l}^{c}f(\boldsymbol{c}_{\nu,l})$ is convex in $\boldsymbol{C}$ as sums of convex functions are also convex. The functions in the constraints given in (20) are linear in $\boldsymbol{C}$ and hence convex. Therefore to find the $\boldsymbol{C}$ which maximise $\mathcal{L}(\boldsymbol{\theta},\boldsymbol{\theta}^{c},\delta)$ (16) we can use the same numerical approach as used to minimise (3). The expectation and maximisation steps are applied iteratively until the the change in value of the Lagrangian function from one step to the next is less than a prescribed threshold. We used a value of 1e-4 for the results in this paper.

**Initialisation

**Depending on the starting conditions the EM algorithm may converge to a local optimum. To overcome this issue we use a form of the random restart approach commonly used with the EM algorithm for Gaussian mixtures (Biernacki et al., 2003) described as follows. To set initial estimates for $\boldsymbol{C}$, each observation $\boldsymbol{\chi}_{l}$ is randomly assigned a class label $g_{l}^{\prime}\in{1,2,\ldots,k}$. Rate function estimates for the ‘random’ classes are obtained using the classification procedure outlined in Section 3 and the basis function estimates which correspond to these are used for the initial estimates for $\boldsymbol{C}$. All mixing weights in $\boldsymbol{\tau}$ are initially set to 1/$k$. The EM algorithm is then initialised from a number of different sets of initial conditions corresponding to different ‘random’ class assignments and the solution with the largest log likelihood is retained. Using this procedure, for all the examples considered in this paper, 3 random restarts were typically plenty to avoid local maxima.

5 Empirical results

In this section we present the results of our methods applied to various synthetic and real world data sets.

5.1 Synthetic data

We first applied our methods to synthetically generated event data with prescribed rate functions. None of the data generating rate functions belong to the class of smooth functions we assumed in our methods. We generated three difference data sets based on the following four rate functions.

[TABLE]

[TABLE]

[TABLE]

For synthetic data set 1 we generated event times for observations with class label 1 (resp. 2) separately using rate $\lambda_{1}(t)$ ($\lambda_{2}(t)$). As $\lambda_{1}(t)$ and $\lambda_{2}(t)$ are both periodic functions this experiment mimics real world examples of point processes which exhibit periodic patterns (e.g. time of stock trades (Engle, 2000) and the times at which emails are sent (Malmgren et al., 2008)).

For synthetic data set 2, $\lambda_{3}(t)$ (resp. $\lambda_{4}(t)$) was used to generate class 1 (2) observations. These rate functions mimic computer network traffic for which a NHPP model with rate function defined by a step function is a suitable model (Paxson and Floyd, 1995).

Finally for data set 3 we considered a four class problem where rate functions $\lambda_{1}(t),\lambda_{2}(t),\lambda_{3}(t)$ and $\lambda_{4}(t)$ were used to generate data in classes 1,2,3 and 4 respectively. The thinning algorithm (Lewis and Shedler, 1979) was used to generate all data. Synthetic event time data for 20 observations for each rate function is shown in Figure 1.

5.1.1 Classification

For each class, event time data for 20 observations was generated. 10 of these formed the training set and the remainder formed the unseen test set. Using the classification procedure outlined in Section 3 we were able to recover the prescribed rate function for each data set (see Figure 1). Furthermore, when each observation in the test data was assigned to the class for which its membership probability (as defined in equation (10)) is maximal, the classification accuracy was 100%.

5.1.2 Clustering

Again, data for 20 observations was generated. Using the clustering method detailed in Section 4 we were able to recover the NHPP rate functions for all data sets (see Figure 1). For each data set, we were able to correctly assign each observation to the NHPP model from which it’s event time data was generated by assigning it to the model for which the membership probability (equation (15)) was maximal.

5.2 Retail transaction data results

Next we introduce our first real data set which is made up of the till transaction times from the 6${}^{\mbox{th}}$ of September 2011 to the 28${}^{\mbox{th}}$ of February 2015 for 74 UK stores (10,000 transactions per store) belonging to a large UK retailer. Half of these stores have been categorised by the retailer as ‘High street’ (HS) which are large stores found in city centres and the remainder as ‘travel’ stores which are situated in airports and railway stations. These categorisations serve as the ground truth for our analysis.

5.2.1 Classification

We used 5 fold cross validation (CV) (Efron and Gong, 1983) to obtain estimates for the classification accuracy, high street store and travel store true positive rates (TPR). We repeated the CV procedure 100 times with different random fold partitions to give a total of 500 numerical experiments. The average classification accuracy, high street store true positive and travel store true positive rates were 0.847, 0.994 and 0.701 respectively. The estimates for the rate functions of each class for one CV run are shown in Figure 2. Both rate functions are periodic with the high street rate function peaking around Christmas which corresponds to the busiest shopping period on the high street. In contrast, the rate function for the travel stores (which are located in airports and railway stations) peaks around summer which coincides with the peak holiday period when holiday makers are travelling via transportation hubs. For our results, the high street store true positive rate was higher than the travel store true positive rate suggesting the classifier has more of a tendency to misclassify travel stores rather than high street stores. Further investigation reveals that the misclassified travel stores are situated in central railways stations in large UK cities including London, Glasgow, Leeds and Liverpool. A possible explanation for the misclassification of these stores is that, although classed as travel stores, the store location means that these stores predominantly serve high street shoppers and thus have a demand profile which is similar to high street stores.

5.2.2 Clustering

Setting the number of models in the mixture to 2, and assigning each observation to the model for which its membership probability is maximal, the average clustering accuracy, high street store true positive and travel store true positive rates were 0.764, 1 and 0.528 respectively. The high street store true positive rate (resp. travel store true positive rate) of the clustering results was higher (lower) than the corresponding classification results. Further investigation in to the cause of this discrepancy revealed that most of the stores in the ‘travel’ cluster are located in airports. Stores categorised by the retailer as travel stores, but with temporal transaction characteristics more similar to high street stores (such as the stores situated in central railway stations in large UK cities mentioned in the previous section) cluster into the ‘high street’ group. The rate functions of the two NHPP models are shown in Figure 2 and clearly the blue rate function peaks around summer and corresponds to travel stores,while the red rate function peaks around Christmas and corresponds to high street stores.

Rate function estimates obtained by setting $k=3$ are also shown in Figure 2. Two clusters, with rate functions labelled ‘airport 1’ and ‘airport 2’ in the plot are comprised almost exclusively of airport stores and have very similar rate functions peaking around the summer. A third cluster, labelled ‘high street and railway’ is made up mainly of high street stores and railway stores located in stations in city centres and its rate function peaks around Christmas. The fact that setting $k=3$ does not lead to more meaningful store partitions, compared to setting $k=2$, supports the notion that the stores fall naturally in to one of two clusters, the first made up predominantly of airport stores and the second of stores located centrally in cities either on the high street or in city centre railway stations.

5.3 Hubway bike share data set

This data set contains bike trip information for the Hubway bike share program located in Boston, Massachusetts and its environs. The data, which covers the period from the 28${}^{\text{th}}$ of July 2011 to the 30${}^{\text{th}}$ of November 2013, contains trip duration time, date and time of the trip as well the the id of the start and finishing bike share stations of cycle journeys. Additionally the user type (casual - ‘24-hour or 72-hour pass user’ or member - ‘annual or monthly member’) was also recorded. In our analysis we ranked each station according to the proportion of station users classed as either casual or members and classified the $\kappa$ stations with the highest casual-member proportions as ‘casual’ stations and the $\kappa$ stations with the lowest casual-member proportions as ‘member’ stations. We only considered stations with more than 5,000 recorded trips to avoid any spurious results caused by stations with only a handful of recorded trips. This left a total of 89 station observations. The times of day at which the bike trips began at each station served as the sample of event times for each station observation.

5.3.1 Classification

We used the same CV procedure as before to obtain average classification accuracy, casual station and members station true positive rates. Raster plots showing the time of day at which bike trips started for both station types and plots showing the average performance of our classification procedure over all test folds as a function of $\kappa$ are shown in Figure 3. The results show that our classification procedure did a very good job of classifying the station observations (the average classification accuracy for $\kappa=20$ is 0.775 for example). We also include a plot showing the estimated rate functions for the casual and member station classes for $\kappa=20$ in Figure 3. The member stations rate function (red curve) has two distinct peaks around 8am and 5pm. These peaks are far less pronounced for the casual stations (blue curve) which has a larger intensity in the middle of the day. The most plausible explanation for these features is that member stations are predominantly used by commuters who cycle to and from work and, as they use the bike scheme every day, are more likely to sign up for membership. Whilst casual stations are more likely to be used by tourists who have the free time to cycle during the middle of the day but, as they use the bike share scheme infrequently, are less likely to become members.

Although the classification results are very good, for the 2 class problem, our classifier had a tendency to misclassify member stations more than casual stations. In Figure 3 is shown the estimated rate functions for the two ‘member stations’ which were misclassified the most number of times during the CV experiments. Rather than two distinct peaks around 8am and 5pm, station 9 has a single peak around 8am and station 15 a single peak around 5pm. This suggests that there may be two distinct sub-types of member stations; one type used mostly in the mornings (most likely located in the suburbs close to the homes of commuters) and another close to the commuters places of work (most likely in the city centre) which are mostly used after 5pm when the commuters return home from work. This motivated us to consider a three class problem consisting of ‘casual stations’, ‘morning member stations’ (member stations where the majority of trips took place before 5pm) and ‘afternoon member stations’ (member stations where most trips took place after 5pm). The results of this experiment are also shown in Figure 3 where it can be seen that the average classification accuracy for the 3 class problem is higher than for the 2 class problem (for $\kappa=20$, the 3 class problem average classification accuracy is 0.829, whilst it is 0.775 for the 2 class problem for example). Furthermore, the rate function estimates for the 3 class problem, which are also shown in Figure 3, capture the different usage patterns of the 3 different station types for the 3 class problem. There is a clear peak in activity for the morning (resp. afternoon) member station class around 8am (5pm), whilst the estimated rate function for the casual station type has more of a uniform profile.

5.3.2 Clustering

The clustering results for the 2 and 3 class problems are shown in Figure 3 and show that the stations classes defined in the previous subsection emerge naturally from the data. For $k=2$, the estimated rate function for one station type has two distinct peaks at 8am and 5pm. This is the ‘member station’ class. The other estimated rate function has a greater intensity during the middle of the day. This is the ‘casual station’ type, although for the clustering results, it has a noticeably higher peak around 5pm when compared to the rate function of the casual station class obtained via the classification procedure for $k=2$. For $k=3$, the casual station class persists and the two morning and afternoon member station sub classes, which were defined in the previous subsection, emerge naturally from the data. Assigning each observation to the class for which its membership probability is maximal, the clustering accuracy and true positive rates were very good (for $\kappa=20$, the clustering accuracy is 0.750 for $k=2$ and 0.889 for $k=3$) and broadly reflect the results seen for the classification task. For $k=2$, the casual station true positive rate exceeds that of the member station true positive rate. However, setting $k=3$ goes some way to remedying this discrepancy and also results in a higher overall clustering accuracy, as it did or the classification results. This indicates a choice of 3, rather than 2, clusters is a more appropriate choice for the clustering problem.

5.4 A&E data set

This data set contains the time of day (00:00:00 to 23:59:59) in every month from April 2011 to the end of December 2014 at which patients arrived at Accident and Emergency (A&E) departments in hospitals throughout England. Each patient was diagnosed as either suffering from poisoning (including alcohol poisoning) or a cardiac condition and these are used as the ground truths for the class labels. For each label there are 45 observations corresponding to each of the 45 months spanning April 2011 to the end of December 2014. Each monthly observation contains 10,000 A&E admission events.

5.4.1 Classification

Using the same CV procedure as outlined in Section 5.2.1, we obtained an average classification and true positive rates for both classes of 1. The estimated rate functions are shown in Figure 4 from which it can be seen that admissions for poisoning peak around midnight. This is most likely due to instances of alcohol poisoning amongst revellers in the late evenings. Cardiac conditions peak in the morning, which is consistent with evidence in the medical literature (Peters et al., 1989, Muller et al., 1989).

5.4.2 Clustering

Estimated rate functions from the cluster analysis are also shown in Figure 4 and these are very similar to the estimates obtained from the classification procedure. In particular, instances of poisoning peak around night time and cardiac conditions reach a high point in the morning. Assigning poisoning and cardiac months to the model for which their membership probabilities are maximal gives a clustering accuracy of 100%.

We have made Matlab code for the classification and clustering procedures outlined in Sections 3-4 and for the synthetic data set results in Section 5.1 available at https://github.com/duncan-barrack/NHPP. The Hubway bike share data set is available at https://www.thehubway.com/system-data. The store data and NHS A&E data are not publicly available.

6 Discussion

In this paper, we have detailed principled approaches for the classification and clustering of observations of event data using NHPP models. Results on synthetic and real data were presented which show the effectiveness of our methods. The focus of this work has been on temporal point process data observed over a fixed interval $\{t:0<t\leq T\}$. Another possibility, particularly suited for periodic data, would be to assume the domain of the data is periodic (e.g. a circle). The approach we have described here should generalise easily to the circular case. Furthermore, our method is not restricted to temporal data and, in principle, can be generalised to multi-dimensional spatial and spatial-temporal point process data.

For the clustering procedure choosing the appropriate number of models in the NHPP mixture may not always be straightforward. In principle, the Bayesian information criterion (Schwarz et al., 1978) or Akaike information criterion (Akaike, 1974) could be used for this. Future work could investigate the suitability of these methods for choosing the number of NHPP components for our clustering procedure.

Acknowledgements

We would like to thank The UK based retailer and NHS England for their support and for providing the data used in our analysis. This work was funded by EPSRC grant EP/G065802/1 - Horizon: Digital Economy Hub at the University of Nottingham and EPSRC grant EP/L021080/1 - Neo-demographics: Opening Developing World Markets by Using Personal Data and Collaboration.