Z-Transforms and its Inference on Partially Observable Point Processes

1 Introduction

Background and issues. Point processes are statistical models for modelling random times at which events of interest take place. In a variety of applications in machine learning, maximum likelihood estimation can be difficult for point process fitting, depending on the specificities of the point process itself. As the likelihood may not be convex in the parameters, a maximization routine would likely converge to merely a local maximum as opposed to the intended global maximum.

One standard strategy used to identify the global maximum involves using several sets of different initial values for the maximization routine. We remark that this plan of action does not solve the problem in its entirety, and a local maximum may still be identified erroneously as the global maximum. This strategy can be taken further: more than one different maximization routines, with very distinct optimization techniques, may be used collectively with several sets of initial values to infer the parameters of interest. If these different maximization routines identify similar points as the potential global maximum, more confidence can be placed in judgements that the identified points are indeed the actual global maximum.

Another possible approach that is frequently used in the literature is the moment-based estimation wherein parameters are obtained by minimizing a measure of the discrepancy between empirical and theoretical second-moment properties. These raw and uncentered moments are usually easy to compute and they reveal important aspects of a distribution. For example, the first four moments tell us something about the information regarding the mean, variance and skewness as well as kurtosis, respectively. Using this information, one can immediately place constraints according to our theory on the so-called location, scale or shape and the tail behaviour of the distribution without specifying a full model.

The same line of thought applies to higher order moments. One might ask whether it is possible to match a larger set of moments, or the entire spectrum of moments (rather than truncating at the first four moments, say) in order to get an even closer representation of the model. If so, what would be a suitable quantity to match? This paper seeks to address these questions and shows that inference is possible through matching a transform of the raw and uncentered moments.

Purpose of the present study. We formulate an alternative estimation method using the Z-transforms for a class of nonhomogeneous point processes whose intensity is driven by a general neural network. The Z-transform is an integral transform in which this transformation 'compresses' the point process into a single function of an argument, η, which we call the dummy variable. This dummy variable η serves as placeholders for the probabilities that determine the distribution of the point process.

Z-transforms are in general an alternative specification of stochastic processes that somehow encode the properties of the distributions into a form that is more convenient for certain kinds of probability calculation. We note that the Z-transform with dummy η is also known in combinatorics as the probability generating function.

This quantity fully characterizes the distribution of the point process in the sense that there exists a one-to-one correspondence between the point process structure and the Z-transform. Given the form of our neural network intensity driving the point process, we exploit the structure of the Z-transform to estimate the parameters of interest. The inference strategy is to match the theoretical Z-transform with its empirical counterpart from the observed samples.

In general, we may arrive at an over-determined situation: that is, one in which an abundance of Z-transform conditions are available. Similarly, there might be too few conditions than parameters to estimate wherein there is insufficient information and the model is under-determined. Through our formulation, we match the entire spectrum of the Z-transform by regulating η so that there is a surfeit of moment equations outweighing the number of parameters as our ideal situation.

Event censoring. These Z-transforms are brought only to illustrate one rather important application in point processes - censoring. Censoring is a phenomenon in which the value of an observation is only partially known. In the present context, the times at which events happen are censored. Some notable examples include medical trials: the subject may be observed only periodically, rather than assessed in a continuous fashion. The medical practitioners would report the number of heartburn per day rather than the exact times at which the heartburn occurs within that day for patients.

Contributions. Our noteworthy contributions are as follows:

1. We present a method of inference via the Z-transform for a specific class of non-homogeneous point processes whose intensity is driven by a general neural network. The idea is to match the theoretical Z-transform in its entirety with its empirical counterpart from the observed samples. This is done by regulating the number of dummy parameters η so that we have as many equations as there are moment conditions. This would ensure a single solution for the parameters.
2. We use this perspective to solve inference problems motivated by events that are censored. Formally, censored data arises when the event of interest cannot be directly observed but is only known to have occurred during an a period of time. In such cases, the only information we have is the number of events over a given period, but the times at which the event happens are not known. Through Z-transforms, we show how to carry out inference in such situations.
3. We explicate a classical connection where the k-th moment for our model exhibit a closed-form solution in terms of Stirling numbers and the Bell polynomial. The predictive formula k = 1 and expressions for higher moments to calculate variances, skewness and kurtosis are readily available.
4. We present experiments in a controlled setting and show that our inference algorithm are able to indeed recover the optimal parameters.

2 Notational Conventions and Prerequisites

Our framework requires some background on temporal point processes. Let us cover this topic briefly.

2.1 Temporal Point Processes

Temporal point processes model the times at which events of interest occur via a stochastic process (N(t))t≥0, where N(t) -N(s) measures the number of events that occur in the time interval (s, t]. We focus on such a process.

Fix some locally integrable λ : R+ → R+. We call (N(t))t≥0 an non-homogeneous Poisson process with intensity λ if

1. N0 = 0 almost surely
2. for any s < t, N(t) -N(s) ∼ Poisson (∫ t s λ(x) dx)
3. for any s < t ≤ s′ < t′ , N(t)-N(s) ⊥⊥ N(t′ )-N(s′ ).

One may view non-homogeneous point processes as the following generative model for event times: given some terminal time T , the number of events N is drawn from a Poisson random variable with mean ∫ t s λ(x) dx, and the N event times are then drawn i.i.d. from a distribution with probability density function given by

3 The Model - Neural Network Intensity

The nature of event counts in many applications is likely to be time-varying. A cursory visual inspection of cyber breaches [Privacy Rights Clearinghouse, 2018] over time (see Figure 1) reveals several periodicities. We remark that there is cyclical variability in the sense that there are some cyclical trends or seasonal trends (i.e., some repeated patterns of cyber breaches).

To this end, we formulate an intensity that is general enough to encapsulate these time varying characteristics, while presenting an inference mechanism via matching the Z-transform to learn the appropriate activation functions. We present a flexible family of a neural network

where τ > 0 being the background intensity driving the point process, g(·), f(·) are non-negative activation functions, b(t, δ) is a decreasing function in t with δ as the parameter, and R is a fixed number which we call the hidden units. This model allows us to control the trend behavior in the intensity function, as can be seen in Figure 2. The parameters to estimate are θ = {τ, δ, (ar, cr, dr)Rr=1 }. Each term f(cr ·t+dr) may be seen as a learned hidden representation of the input time. By choosing R to be sufficiently large, we can model complex non-linear hidden periodicities within the rate of arrivals governing the number of events. A standing assumption that we make is the follows:

This intensity capture a wide range of other well known intensities. Some examples include

1. Polynomials.

with the output activation function taking the identity function.

2. The homogeneous Poisson process.

with the output activation function g being the identity function.

3. The p-th oscillation intensity.

with the decaying function ar = 1, b(t, δ) = e-δt . Figure (2) illustrates the dampening of this intensity function by controlling the parameter δ > 0 and the p parameter.

4. Sigmoid intensity.

with the output activation function taking the identity function.

Existing work. The idea of using neural network intensities for general point processes is not new; several authors [Kim et al., 2011; Du et al., 2015; Choi et al., 2016; Mei and Eisner, 2016; Du et al., 2016; Xiao et al., 2017; Wang et al., 2017] have recently proposed the use of recurrent neural networks in conjunction with point processes with random intensities. These more complex models require Bayesian inference and likelihood methods, which prevents scalability. Our intensity function is inspired by, and a generalization of, the one proposed in [Menon and Lee, 2017].

3.1 The Z-Transform

Here we compute explicitly as possible the Z-transform for our model. The Z-transform is defined to be

for η ∈ (0, 1). As mentioned earlier, this is also known as the probability generating functional [Cox and Isham, 1980] or at times known as exponential tilting. We remark that this quantity uniquely determines the distribution of the entire point process [Kallenberg, 2006].

Proposition 1 Let N be a point process and λ(·, θ) be its associated intensity process. Let η ∈ (0, 1). Then

provided the assumption in equation (3) holds true.

Proof. The proof follows similar lines to that of Watanabe's Charaterization Theorem, confer [Br´emaud, 1981]. For a fixed T < ∞, the integral in equation (1) is finite due to standing assumption as stated in equation (3). Note however that the original Watanabe's Characterization Theorem deals with expression of the form

but evaluating the quantity of interest in equation (8) requires one to only perform an exponential tilting tweak of a quantity, so the extra details of the proof is elided.

3.2 The First k Moments

The flexibility of theoretical Z-transform is borne out in evaluating the first k moments. As pointed out in [Pitman, 1997], one can deduce the following

where the { kz } are known as the Stirling numbers of the second kind that counts the number of ways to partition a set of k elements into z sets and B is the Bell polynomial. This is an interesting connection because, one might be interested in finding the higher order moments of our model, after parameter estimation of ˆ θ and this can be evaluated in Python by using the mpmath.bell() function.

The predictive formula. The predictive formula for a period of s′ , t′ would take the following form:

This predictive formula gives the expected number of events following the learned parameters ˆ θ using the procedures of inference deferred to Section 4.

4 Estimation Procedure

We spell out an optimization technique to learn the parameters θ of our point process driven by a general neural network intensity function. Recall that our intention is to match a spectrum of theoretical moments (not just the mean, variance, skewness or kurtosis) through the analytic Z-transform of a neural point process as computed in Proposition 1.

It is natural to assume that the following holds: if θ0 is the true value of the parameter θ, then the theoretical Z-transform moments and the empirical moments must match.

Let ∆j(t, s) be the observed number of events in sample j during the observation period (s, t]. Then the empirical Z-transform is given by φm(η, t, s) where

From Proposition 1, the theoretical Z-transform (repeated here for convenience) for our neural point process model is

The parameters are obtained by minimizing the relative entropy between the Z-transform and its empirical counterpart, i.e.

where p denotes the normalized empirical probability and q the normalized theoretical probability of the Z-transform. We need the normalization since we have to convert the values of Z-transforms and its empirical counterpart to probability distributions.

Sigmoid intensity λσ(·, θ). We illustrate this with a concrete example. Suppose we fix our neural point process model as having a sigmoidal flavored intensity, i.e.

The function has the property that the initial stage of intensity growth is approximately exponential; then, as saturation begins, the value slows, and at maturity, it terminates. Observed that there are three parameters to be estimated, namely {τ, c, d}. The theoretical Z-transform associated with the sigmoid intensity λσ(·) is calculated to be

Note in this particular example, there are three unknowns. The normalized distributions of q and p are calculated as follows:

Algorithm 1 Inference of parameters via Z-transforms

Input: Fix time points (s, t), m observed number of events for fixed times and ∆j(s, t) for j = 1, 2, . . .m. Define the intensity specification λ(t).

• 1: Given the observed number of events from s to t, ∆j(s, t), of m number of paths, compute the empirical moments:

• 2: (a) For the intensity specification that leads to a closed form of Z-transforms, compute this quantity via:

• (b) For the intensity specification that does not inherit a closed form Z-transform, compute the theoretical moments approximately with an aid of a using symbolic programming and numerical integration.
• 3: Compute the normalized distributions q and p via expressions in equations (18) and (19) respectively.
• 4: Minimize the relative entropy or mean squared error between the theoretical and empirical Z-transforms as in equation (15) or equation (21) using standard optimization packages whilst enforcing strict constraints such that τ > 0, δ > 0 and g(·), f(·) being non-negative activation functions.

Hence, the estimated parameters associated with λσ is then constructed by minimizing the relative entropy between the theorerical Z-moment and its empirical version calculated from the observed samples, i.e.

Another way to estimate λσ is to minimize the mean squared loss between the theoretical Z-transform and the empirical version as follows:

Note that this is the standard loss used in the GMM framework [Hall, 2005].

Overdetermined and underdetermined systems. Let the dimension of θ be p and r be the number of conditions. In general there might be more moments conditions than parameters to estimate r > p, i.e., the model is overidentified. If r < q, then there is insufficient information and the model is underidentified.

The role of η in the Z-transform. Rather than matching the first k moments viz. mean, variance, skewness and kurtosis, say, we match the entire structure of the Z-transform by regulating the number of different values for η until there is a surfeit of moment conditions compared to number of unknowns.

Connections to literature. In most of the point process literature, the moment matching is carried out by matching the first few moments as well as the tilted moments including the autocorrelation functions, [Zhang et al., 2012; Yu, 2004; Gerhardt and Nelson, 2009; Heath et al., 2013; Diggle et al., 2013; A¨ıt-Sahalia et al., 2015; Da Fonseca and Zaatour, 2014]. Similar techniques were extended to inferring non-parametric forms of point processes that relies on second-order and third-order integrated cumulants [Achab et al., 2017].

Depending on the form of the objective and the iterative specificities, these methods can be seen as falling into the framework of the generic method of Generalized Method of Moments for estimating parameters in probabilistic models [Hall, 2005] and are used in other applications of machine learning [Anandkumar et al., 2012; Soufiani et al., 2013].

5 Experiments

5.1 Parameter Recovery on Fictitious Data

Our inference algorithm is tested on synthetic data generated from our model. We validate our theoretical analyses by illustrating the feasibility of using the Z-transform matching technique for our neural point process model.

Preliminary setup. For fixed time points (s, t), we use the modified thinning method, confer [Ogata, 1981]; to generate m paths from a neural point process with known parameters θ∗ . For each of this sample j, where j = 1, 2, . . . ,m, we record the total number of events and we compute the parameter estimate ˆ θ via minimizing the procedures set forth in Section 4.

Intensity functions. We repeat this for 5,000 independent samples from the point process for two distinct intensity functions, i.e., 1. linear intensity which takes the form

and 2. the oscillating intensity, whose intensity diminishes as time goes by taking the following form:

Note that this is a special case of that given in equation (6).

We set the parameters θ∗ for the linear intensity as θ∗ 1 = (τ, a) = (0.5, 0.85) and the oscillating intensity as θ∗ 2 = (τ, δ) = (3.0, 0.15).

Results. Tables 1 and 2 confirm that the Z-transform method have commensurate accuracy to the ground truth θ∗ 1 and θ∗ 2 . This reassures that when the Z-transform method is used, these inference procedures operate sensibly, and indeed recover the optimal parameters asymptotically.

With regards to the linear intensity (see Table 1), this intensity may seem trivial but this is an example where both the MLE as well as the our Z-transform method have corresponding accuracy to the ground truth θ∗ . However for the oscillation intensity (see Table 2), we see that the MLE method turns out to be problematic, possibly owing to the non-convexity and periodic nature of intensity function. We report that the calibrated parameters by matching the first two moments for linear and oscillation intensities are (τ, a) = (0.4000, 0.8997) and (τ, δ) = (0.0458, 0.6996), respectively. Also in addition, note that the oscillation intensity function has the form constant+exp(-δt) · cos(2π t) where the intensity decays off as time t ↗∞. In this particular situation, the effect of the constant dominates the generation of the events. This could be the reason the process of recovering the ground truth for δ is a little challenging, as evident from Table 2.

(a) Calgary

Table 3: Comparison of compute time (in seconds).

As a caveat on the applicability of the our method, we remark that the expression of exponential in equation (9) may at times cause numerical issues since for very large numbers of events, the resulting value for the transform is unrepresentably tiny.

5.2 Event Counts Predictions on Real Data

We compare the various methods on two real-world datasets: breaches data, henchforth known as Breaches, comprising times of security violation in a period of 2005 to 2018 [Privacy Rights Clearinghouse, 2018]; and Calgary, comprising of access log to the University of Calgary from October 1994 to October 1995 [The Internet Traffic Archive, 2018]. Previous studies have utilized standard time series models to fit and predict potential cyber attacks [Kumar et al., 2006]. For each dataset, we aim to predict the number of events happening in some specified future time window.

Preliminary setup. For each of these datasets, we split the recorded number of events into a train and test set.We then carry out predictions of event counts over periods with different granularities. The precise split methodology varies for each dataset are as follows: For the Breaches dataset, we select the first 3,926 data points as training set to train the models and predict the next 1 to 4 data points (in days). For Calgary, the first 5,000 data points are selected as training set and we predict for the next 20 minutes. The prediction is done via these intervals using the Bell function formula explicated in Section 3.2.

Results. Figure 3 confirms that on real-world data, the performance of our model seemed competitive compared to standard baseline models: AR(1), GARCH(1, 1), as well as the Gaussian Process Regression (GPR). In Table 3, we show the training times of our method compared to AR(1), GARCH(1, 1), and GPR. The results seemed to indicate that there are benefits over existing models in performing short term predictions. The arguments as to why these might be the case can be found in the recent work of [Menon and Lee, 2017]. Our method is the fastest, while GPR is the slowest. All other methods including our method work well on big datasets; the GPR method seems unable to fit a training set with more than 6,000 points, as it faces a memory usage issue (using the scikit-learn implementation in a Core i7 machine with 8GB RAM); future work may compare against recent sparse Gaussian Process models or online approaches for time-series (e.g., [Soh and Demiris, 2015]).

6 Conclusion

We proposed a method of drawing inferences via the Z-transform for a point process whose arrival rate takes a neural network. The idea is to match the theoretical Z-transforms with its empirical counterpart of the observed samples. We further used this perspective to investigate the inference mechanism in situations wherein events are censored, i.e., the only information we have is the number of events over a given period, but the times at which the event happens are not known. In future work, exploring similar ideas for larger class of point processes whose governing intensity functions are in fact random (confer e.g. [Linderman and Adams, 2014; Flaxman et al., 2017; Ding et al., 2017]) would be of interest. Another area we expect to make progress in is the study that concerns properties of Z-transform estimators.

Acknowledgments

Young Lee wishes to thank Aditya K. Menon for his careful and critical reading of our manuscript.

This research is supported by the National Research Foundation, Prime Minister's Office, Singapore under its National Cybersecurity R&D Programme (NRF2015NCR-NCR002-001).