[2]Hoai Nguyen Xuan These authors contributed equally to this work.

1]Department of fundamental science, Academy of Cryptography Techniques,
No.141 Chien Thang, Hanoi, 10000, Vietnam

[2]Faculty of Information Technology, HUTECH University, No. 475A Dien Bien
Phu, Ho Chi Minh City, 72308, Vietnam

# Learning Generalization and Regularization of Nonhomogeneous Temporal
Poisson Processes

Son Nguyen Van<EMAIL_ADDRESS><EMAIL_ADDRESS>[  *

###### Abstract

The Poisson process, especially the nonhomogeneous Poisson process (NHPP), is
an essentially important counting process with numerous real-world
applications. Up to date, almost all works in the literature have been on the
estimation of NHPPs with infinite data using non-data driven binning methods.
In this paper, we formulate the problem of estimation of NHPPs from finite and
limited data as a learning generalization problem. We mathematically show that
while binning methods are essential for the estimation of NHPPs, they pose a
threat of overfitting when the amount of data is limited. We propose a
framework for regularized learning of NHPPs with two new adaptive and data
driven binning methods that help to remove the ad-hoc tuning of binning
parameters. Our methods are experimentally tested on synthetic and real-world
datasets and the results show their effectiveness.

###### keywords:

Nonhomogeneous Poisson Process, Learning Generalization, Arrival Rate
Regression

##  1 Introduction

Stochastic processes allow us to model uncertain quantities that evolve over
time. One of the most fundamental and important of such processes is the
Poisson process [1, 2, 3, 4] as it is appropriate for modelling a large number
of completely random phenomena. Therefore, the Poisson process has been widely
used in various real-world applications such as taxi routing [5], queueing
systems [6], and e-mail arrival counting [7, 4].

Because of the importance and popularity of Poisson processes, estimation of
them from data has been extensively studied in the field of statistics and
operations research, especially for the one-dimension case (called Temporal
Poisson Process). While estimation of homogeneous Poisson processes (HPP),
where the rate is constant, is trivial, the problem of estimating
nonhomogeneous Poisson processes (NHPP), where the rate is a function of time,
is considered more challenging [6, 4, 8]. Moreover, all of the estimation
methods in the literature so far have relied on dividing observed data into
bins (binning method), very often with equal length [3]. While it is
reasonable to use binning method for statistical estimation with infinite
amount of data as it helps to reduce the estimation difficulty by transforming
a NHPP into a HPP (see further explanations in section 4), it could create
some serious problems when working with finite amount of data (which is
usually the case for real-world applications). Firstly, when the amount of
data is finite, dividing into bins might not guarantee that the data in each
bin still follows assumed distributional assumptions, in this case - Poisson
property. For instance, in all the works on estimation of NHPPs in the
literature that we have studied, none of them even tested if the observed data
used in their experiments actually came from a Poisson process (except for
synthetic datasets coming from simulating some predefined hypothetical Poisson
processes). Some even explicitly acknowledged this concerned research hole
[4]. Secondly, with finite amount of data, the estimation becomes a learning
problem, where the central issue is generalization, i.e how well the
estimated/learned model could predict the future/unseen data. To the best of
our knowledge, we have not seen in the literature any work that systematically
formulate and solve the problem of estimation of NHPPs from finite amount of
data as a learning generalization problem. Thirdly, how the binning should be
done so that it could help to solve the first and second aforementioned
problems.

In this paper, we formulate the problem of estimating temporal NHPPs from
finite amount of data as a learning generalization problem. Our main
contributions are:

  * •

We mathematically prove that for estimation/learning NHPPs from finite amount
of data, data binning could increase the probability of overfitting (i.e the
learned model is good on observed/training data but bad on future/ unseen
data).

  * •

We propose two adaptive and data-driven binning methods that help to
regularize the learning of NHPPs as well as automatically obtains the number
of bins so as to remove the ad-hoc tuning of this parameter for the
estimation/learning process.

The remaining part of this paper is organized as follows. The related work is
briefly given in Section 2. Section 3 provides some backgrounds on
mathematical concepts and results used later in the paper. The problem
description is shown in Section 4. Our NHPPs learning framework and the two
novel binning methods for learning regularization of NHPPs are detailed in
Section 5. We present and discuss our experimental results in Section 6. The
paper is concluded in Section 7.

##  2 Literature Review

There are a number of approaches to estimate NHPPs via estimating their rate
functions in nonparemtric and parametric manner. For nonparametric approaches,
[9, 10, 8] proposed simple estimators for estimating the cumulative rate
function of NHPPs on the time interval from historical data. Their methods do
not require the assumption of any parameters or weights for the rate function.
The authors in [9] considered interval estimators of the cumulative rate
function. The expected number of events for each interval is a piecewise-
constant rate function. In [10], the authors extended the work of [11] to
solve one or more realizations (data observations) that include overlap ones
on an interval. Recently, the authors in [12] exploited analogies between
point process and time series models. An auto-regressive moving-average model
for NHPP process is proposed to disentangle cluster mechanisms in continuously
observed count data of an NHPP. Their algorithm overcomes the limitation of
traditional estimation by allowing for a variety of parametric specifications,
which need not be confined to the Markovian case, nor to fixed intensity. An
obvious deficiency of these methods is that their models require every arrival
time of events in every practical realization be stored in memory to allow the
generation of simulated realizations. The works in [8, 2] also have similar
weakness. Other nonparametric estimators of the rate function were introduced
in [13, 1]. In [13], the authors proposed a nonparametric estimator using
wavelets and adapted the nonnegative wavelet estimation of general density
functions proposed in [14] to estimate the rate function. This approach’s main
advantage is that it can be applied to most nonstationary processes that may
or may not show a long-term trend or periodic emergence. In this approach, the
rate function is expressed as a nonnegative linear combination of nonnegative
wavelets. However, a drawback is that they require more constraints on the
rate function than needed to achieve nonnegativity. As a result, there are
several limitations of their procedure in constrained convex nonlinear
programming problems [1].

Most of parametric estimators proposed in the literature have a fixed number
of parameters, so they do not suffer the storage problems when the number of
observed realizations (data) increases. Moreover, when the true rate function
is known to be continuous, piecewise continuous functions could be formed for
fitting the rate functions 111That naturally leads to the binning of the
observed data.. The piecewise function could be linear or nonlinear such as
cubic splines [7, 4]. In these works, the author formulated the problem of
estimation of a NHPP as a semidefinite programming problem (SDP), where the
estimated rate functions are required to be nonnegative. However, they also
doubted if the approach could be practical for large and real world problems.
Their approach do not scale if either the degree or the number of variables
increases, essentially due to intractable numerical calculations. In the
standard SDP formulations, the condition numbers of the feasible dual
solutions are exponential in the degree [15]. Theoretically, every continuous
rate function can be approximated by piecewise continuous functions. However,
when the number of observed data (realizations) is limited, an empirical
estimator might not be asymptotically close to the true rate function if it
does not lie within the assumed parametric class [1].

To our knowledge, almost all the works in the literature have been on the
estimation of NHPPs, where the estimators will converge to the true rate
function when the amount of observed data goes to infinity. None of them has
focused on the case of learning a NHPP from finite data for generalization
over future and unseen data.

##  3 Background

###  3.1 Poisson Distribution

The Poisson distribution, introduced by Siméon Denis Poisson [16],
characterizes random variables $\mathcal{X}$ used for counting the number of
events occurring in an interval. Let the average number of events in an
interval be $\lambda$, $\mathcal{X}$ has Poisson distribution if its possible
values are positive integers and the probability of observing $\xi$ events
($\xi\in\mathbb{N^{+}}$) in the interval is as follows:

$\displaystyle P(\xi\textit{ events in the
interval})=\frac{e^{-\lambda}\lambda^{\xi}}{\xi!}$

For the Poisson distribution, the first two moments, expectation and variance,
equal to $\lambda$ [17]:

$\displaystyle\mathbb{E}[\mathcal{X}]=\Sigma_{\xi\geq
0}\frac{e^{-\lambda}\lambda^{\xi}}{\xi!}\xi=\lambda$
$\displaystyle\textnormal{Var}[\mathcal{X}]=\Sigma_{\xi\geq
0}\frac{e^{-\lambda}\lambda^{\xi}}{\xi!}\xi^{2}-\lambda^{2}=\lambda$

We shall denote $\mathcal{X}$ as $Poisson(\lambda)$.

###  3.2 Temporal Poisson Process

One of the most fundamental and widely-applied stochastic point processes is
the temporal Poisson process used for characterizing discrete event sets in
continuous time. It is usually utilized for modelling the situations where we
want to count the occurrences of certain events that appear entirely at
random, but with a certain rate (without a particular structure) [18].

####  3.2.1 Temporal Homogeneous Poisson Process

Let $\\{\mathcal{N}(t)$, $t\in\mathbb{R}^{+}\\}$ be a counting process for
counting the occurrences of events that occurred from time $0$ up to and
including time $t$. Given $0<t_{1}<t_{2}<\infty$,
$\mathcal{N}(t_{2})-\mathcal{N}(t_{1})$ is the number of events occurred in
the time interval $(t_{1},t_{2}]$.

######  Definition 1.

[19] $\\{\mathcal{N}(t)$, $t\in\mathbb{R}^{+}\\}$ is said to have independent
increments if for any selection of distinct time-points
$0<t_{1}<t_{2}<\dots<t_{m}$, the random vectors
$\mathcal{N}(t_{1}),\mathcal{N}(t_{2})-\mathcal{N}(t_{1}),\dots,\mathcal{N}(t_{m})-\mathcal{N}(t_{m-1})$
are all independent.

######  Definition 2.

[19] $\\{\mathcal{N}(t)$, $t\in\mathbb{R}^{+}\\}$ is said to have stationary
increments if for all $0\leq t_{1}\leq t_{2}<\infty$,
$\mathcal{N}(t_{2})-\mathcal{N}(t_{1})$ has a distribution that only depends
on the length of the time interval $(t_{2}-t_{1})$.

There are several equivalent definitions for a Poisson process. We present
here the simplest one [19].

######  Definition 3.

$\\{\mathcal{N}(t)$, $t\in\mathbb{R}^{+}\\}$ is called the temporal
homogeneous Poisson process with arrival rate $\lambda$ if:

  * •

$\mathcal{N}(0)=0$.

  * •

$\mathcal{N}(t)$ has independent increments.

  * •

$\mathcal{N}(t_{2})-\mathcal{N}(t_{1})\sim Poisson(\lambda(t_{2}-t_{1}))$ for
any two distinct times $0\leq t_{1}<t_{2}<\infty$.

From the above definition, we could see that a homogeneous Poisson process has
both stationary and independent increments.

####  3.2.2 Temporal Nonhomogeneous Poisson Processes

Let $\\{\mathcal{N}(t)$, $t\in\mathbb{R}^{+}\\}$ be a counting process,
$\\{\mathcal{N}(t)$, $t\in\mathbb{R}^{+}\\}$ is a NHPP with arrival rate
$\lambda(t)$ if:

  * •

$\mathcal{N}(0)=0$.

  * •

$\mathcal{N}(t)$ has independent increments.

  * •

$\mathcal{N}(t_{2})-\mathcal{N}(t_{1})\sim
Poisson(\int_{t_{1}}^{t_{2}}\lambda(u)du)$ for any two distinct times $0\leq
t_{1}<t_{2}<\infty$.

We notice that, contrary to the homogeneous case, a NHPP does not have
stationary increments [19, 20].

###  3.3 Distribution Tests

This section gives the descriptions of two main statistical tests for Poisson
distribution. First, the well-known Kolmogorov-Smirnov test (KS) is
introduced. Second, we describe an alternative statistical test based on the
KS test after performing a logarithmic transformation of the data as proposed
in [6].

####  3.3.1 Standard Kolmogorov-Smirnov Test

Suppose that the null and alternative hypotheses of interest are:

$\displaystyle H_{0}:$ The data come from a Poisson distribution.
$\displaystyle H_{1}:$ The data do not follow the Poisson distribution.

The KS test procedure compares the observed data with a specified theoretical
distribution. This test determines the maximum distance between the empirical
distribution function $F^{*}$ of the samples and the theoretical cumulative
distribution function (CDF) of the known distribution $F$, where $F^{*}$ for
$m$ independently and identically distributed (i.i.d.) observations
$\mathcal{X}_{i}$ is defined as:

$\displaystyle
F^{*}(x)=\frac{1}{m}\sum_{i=1}^{m}\mathcal{I}_{(-\infty,x]}(\mathcal{X}_{i})$

where the indicator function $\mathcal{I}_{(-\infty,x]}(\mathcal{X}_{i})$ is
equal to 1 if $\mathcal{X}_{i}\leq x$ and 0, otherwise. The KS statistic for a
given CDF $F(x)$ is calculated as follow:

$\displaystyle D=\sup_{x}|F(x)-F^{*}(x)|$ (1)

The null hypothesis is rejected if $D$ is greater than the critical value
$D_{m,\epsilon}$ obtained from a formula, where $\epsilon$ is the significance
level. There are several variations of this formula in the literature that use
somewhat different scalings for $D$ and critical regions [21, 22, 23], etc. As
in familiar situations, we use the formula
$D_{m,\epsilon}=\sqrt{(-0.5\ln(\epsilon/2)/m}$ in [23].

####  3.3.2 Log Test

For each subinterval $[l,uk]$, we define:

$\displaystyle X_{i}=-(m+1-i)\log(\frac{u-l-t_{i}}{u-l-t_{i-1}}),\textnormal{
}1\leq i\leq m$

where $m$ denotes the total number of events and $t_{i}$ is the arrival time
of the $i^{th}$ event.
Under $H_{0}$, within each given interval, $\\{X_{i}\\}$ will be independent
standard exponential variables. Therefore, any standard test for the
exponential distribution can be applied to test the null hypothesis. For
convenience, we use the KS test. Equation (1) can then be applied using the
exponential CDF $F(x)=1-e^{-x}$. The KS test is widely used for computational
convenience, even though it may not have the most significant possible power
against the alternatives of most interest. Compared to the KS test, [24]
showed the Log Test has higher power, i.e., that it has much greater power
against these alternatives. Therefore, this test is used for our proposed
method in this paper.

###  3.4 Statistical Learning Theory

In this section, we summarize some main concepts and results of the
statistical learning theory given in [25] that will serve as the foundations
for the problem of learning NHPPs.

####  3.4.1 Risk Functional Minimization

Given a set of observations $z_{1},z_{2},\dots,z_{m}$ coming from an unknown
probability distribution $F(z)$, many statistical learning problems based on
these observed data could be formulated as minimizing the risk functional

$\displaystyle R(\alpha)=\int Q(z,\alpha)dF(z)$ (2)

where $Q(z,\alpha)$ is a specific loss function, $\alpha$ is a parameter in
the set of parameters $\Lambda$.

####  3.4.2 The Problem of Regression Estimation/ Learning

One of the common learning problems is regression estimation, which could be
defined as a risk functional minimization problem. Given $G$ as the generator
of vector $x$, where $x\in\mathcal{X}\subset\mathbf{R}^{d}$ is i.i.d.
observations according to some unknown (but fixed) probability distribution
function $P(x)$. The learning machine observes $m$ pairs
$(x_{1},y_{1}),(x_{2},y_{2}),\dots,(x_{m},y_{m})$ where, $(x_{1},\dots,x_{m})$
and $(y_{1},\dots,y_{m})$ are the input vectors and the supervisor’s
responses, respectively. Since the input vectors drawn randomly and
independently from $P(x)$, the supervisor’s responses are obtained at random
from $P(y|x)$. In this case, there exists a joint (but unknown) distribution
function $P(x,y)$. The primary purpose of the learning machine is to construct
an appropriate approximation, which will be used for predicting the
supervisor’s response $y_{i}$ on any input vector $x_{i}$ generated by the
generator $G$.

The approximating function is chosen from a given set of functions
$\overline{F}=\\{f(x,\alpha)\emph{ }|\emph{ }f\in
L_{2}(P),\alpha\in\Lambda\\}$. To choose the best function, we measure the
loss or difference between the supervisor’s response to a given input and the
learning machine’s response. Let the function of conditional expectation is
given by:

$\displaystyle r(x)=\int ydP(y|x)$

This function is called the regression function. We assume that:

$\displaystyle\int y^{2}dP(x,y)<\infty,$ $\displaystyle\int
r^{2}(x)dP(x,y)<\infty$

Since $f(x,\alpha)\in L_{2}(P)$, the minimum of the risk function

$\displaystyle R(\alpha)$ $\displaystyle=\int(y-f(x,\alpha))^{2}dP(x,y)$ (3)

is attained at the regression function $r(x)$ (or the closest to $r(x)$) [25].

####  3.4.3 Empirical Risk Minimization (ERM) Principle

A learning machine is capable of implementing a set of functions
$\overline{F}$. The problem of learning is to choose a function $f(x,\alpha)$
from $\overline{F}$ that minimizes the risk function (3). In order to define
the regression problem, we introduce an $(d+1)$-dimensional variable $z=(x,y)$
and $z_{1},z_{2},\dots,z_{m}$ come from an unknown probability distribution
$P(x,y)$

$\displaystyle z_{i}=(x_{i},y_{i}),\emph{
}x_{i}\in\mathbf{R}^{d},y_{i}\in\mathbf{R},i=1,\dots,m$

We cannot compute $R(\alpha)$ directly but only evaluate the empirical risk

$\displaystyle R_{emp}(\alpha)=\frac{1}{m}\sum_{i=1}^{m}Q(z_{i},\alpha),\emph{
}\alpha\in\Lambda$ (4)

in the situation where the probability distribution $P(x,y)$ is unknown and
the only available information is contained in the set of the observed data.
The set of loss functions $Q(z,\alpha),\alpha\in\Lambda$ is of the form:

$\displaystyle Q(z,\alpha)=(y-f(x,\alpha))^{2}$

The Glivenko-Cantelli theorem [25] provides the theoretical foundation for the
ERM principle in the asymptotic case (i.e. when the number of observations
goes to infinity). It explains why ERM should give good results for learning
problems with a large sample size. However, for a limited amount of empirical
data, in general, a small value of empirical risk does not guarantee a small
value of the actual risk. That leads to the problem of generalization of ERM
based learning machines with small and limited training data [25].

####  3.4.4 Generalization of learning machines

As mentioned in the previous subsection, a small empirical risk (training
error) value does not guarantee a small actual risk (true error) for a
learning machine trained on limited data. It means that a learning machine
could achieve small errors on a (limited size) training dataset but perform
poorly on new (future and unseen) data. In that case, the learning machine is
called overfitting. To avoid overfitting, i.e. to generalize the good
performance on the training data to future and unseen data is the most
desirable property that a learning machine should have. This property is
called learning generalization.

In practice, the generalization capacity of a learning machine is usually
measured by comparing the errors of the learning machine on two independent
datasets - the training dataset (used in the training phase) and the test
dataset (used in the testing phase). In theory, it could be controlled by
using the constructive risk bound given in the statistical learning theory. In
particular, Vapnik [26] proved that with probability at least $1-\eta$,
$(0<\eta<1)$ simultaneously for all functions from the set of totally bounded
functions $0\leq Q(z,\alpha)\leq B,\alpha\in\Lambda$, with the finite VC
dimension $h$, the inequality

$\displaystyle R(\alpha)\leq
R_{emp}(\alpha)+2B\xi\left(1+\sqrt{1+\frac{R_{emp}(\alpha)}{B\xi}}\right)$ (5)

holds true, where

$\displaystyle\xi=\frac{h[\ln{(2m/h)}+1]-\ln{(\eta/4)}}{m}$ (6)

It can be seen from inequality (5) that to ensure that the true risk is small
(i.e. good learning generalization), both empirical risk (training error) of
the learning machine and the VC dimension of the function (hypothesis) space
must be small. In other words, for a guaranteed generalization performance of
a learning machine trained on a limited set of data, the training error should
be small and the learning function should be chosen from a set of functions
with as small VC dimension as possible. However, these two requirements are
contradictory since a hypothesis space as a low VC dimension set of functions
could be too narrow to contain the functions that could achieve small
empirical risk values (in this case, the learning machine is called
underfitting). Therefore, a successful statistical learning scheme should
balance between the choice of the hypothesis space and the chance to get a
function in that hypothesis space with small empirical risk (training error).
One of such schemes is the structural risk minimization principle.

####  3.4.5 Structural risk minimization principle

The structural risk minimization (SRM) principle, proposed in [25], is an
inductive principle for model selection used for learning from finite training
datasets. It describes a general way of capacity control that provides a
trade-off between the chosen model’s complexity and the empirical error.
Let us define a structure of nested subsets:

$\displaystyle\mathcal{H}_{1}\subset\mathcal{H}_{2}\subset\dots\subset\mathcal{H}_{r}\subset\dots$

where, $\mathcal{H}_{k}=\\{Q(z,\alpha_{k})\emph{ }|\emph{
}\alpha_{k}\in\Lambda_{k}\\}$. A structure is called admissible if it ensures
that the capacity $h_{k}$ of $\mathcal{H}_{k}$ is less than $h_{k+1}$ of
$\mathcal{H}_{k+1}$ (see 25 for more details). The SRM principle states that
the subset and the function that minimizes the empirical risk within such
subset (e.g., $\mathcal{H}^{*}$ in Fig. 1), so as to minimize risk bound,
yields the best overall generalization performance. Fig. 1 gives the graphical
explanation of the SRM principle.

Figure 1: Structural risk minimization principle of Vapnik [25]

####  3.4.6 Regularization methods

The regularization method was proposed by A. N. Tikhonov [27]. This method
discourages learning a more complex or flexible model to avoid the risk of
overfitting and increase generalization. Regularization can be as simple as
shrinking or penalizing large coefficients, often called weight decay. The
idea of regularization is to find a solution for the empirical risk functional
(4) that is modified by adding the shrinkage quantity.
Tikhonov regularization minimizes the functional:

$\displaystyle R_{reg}(\alpha)=R_{emp}(\alpha)+\gamma*\textit{Reg}(\alpha),$

where $\gamma$ is called _regularization parameter_ which controls the
importance of the regularization functional $\textit{Reg}(\alpha)$. The
regularization functional $\textit{Reg}(\alpha)$ is chosen to impose a penalty
on the complexity of the model. For the case of Hilbert space, we typically
choose $\mathcal{L}_{2}$-regularization that is differentiable:

$\displaystyle\textit{Reg}(\alpha)=||\alpha||^{2}_{2}$

The coefficients are estimated by minimizing the functional $R_{reg}(\alpha)$.
Therefore, regularization is the process that regularizes or shrinks the
coefficients towards zero. This straightforward solution is one of the main
reasons for the use of Tikhonov regularization in the Hilbert space setting.
Ivanov regularization is yet another method that solves:

$\displaystyle\min_{\alpha}R_{emp}(\alpha),$ $\displaystyle s.t.\emph{
}\alpha\leq\tau$

The Ivanov regularization allows us to handle the two main forces directly: on
the one hand, the depreciation of the functional $R_{emp}$ and, on the other
hand, the control of the hypothesis space, where the hypothesis minimizing the
risk is chosen from.

##  4 Problem Description

For learning a NHPP, it is of primary importance to accurately estimate the
underlying rate function since the entire NHPP could be generated using its
rate function. For instance, the author in [28] specified that realizations of
$\mathcal{N}(t)$ on a specified finite subset $S$ could be generated by
sampling independently on each of the Poisson variates in $\mathbf{R}^{+}$. He
proposed a two-step procedure to generate a NHPP using its rate function
$\lambda(t)$ on a bounded subset $S$. Therefore, for learning a NHPP, our main
focus is to learn its rate function from limited observed data.

###  4.1 Problem of Learning Temporal NHPPs

Given a temporal NHPP $\mathcal{N}(t)$ (as defined in 3.2.2) with a fixed but
unknown rate function $\lambda(t)$, consider time series (the full history of
the NHPP) of a NHPP that are observed
$\mathcal{F}_{[a,b]}=\\{t_{1},t_{2},\dots,t_{m}\\}$ ($a\leq t_{i}\leq
b,\forall i=\overline{1,m}$). The observed dataset $D$ comprises of pairs
$(t_{i},y_{i})$, where $y_{i}$ is the number of events occurred in the time
interval $(t_{i-1},t_{i}]$, i.e., $\mathcal{N}(t_{i})-\mathcal{N}(t_{i-1})$
($\forall i=\overline{2,m}$). As aforementioned, the problem of learning
$\mathcal{N}(t)$ is reduced to learning the rate function $\lambda(t)$ from
the finite dataset $D$, which could be cast to the problem of function
regression as in 3.4.2. Under the setting of the ERM principle, the problem is
to choose a function $\lambda(t,\alpha)$ from a function set $\overline{F}$
that minimizes the empirical risk functional (4).

###  4.2 Generalization of Learning

To approximate the rate function, a popular approach in the literature is
function estimation with infinite data, where a separation between training
and testing phase is not necessary. In an estimation process, a natural first
step is data binning, a frequently-used technique in recent literature [6, 28,
4, 29, 30]. The common estimation process using data binning is given in Fig.
2.

Figure 2: Common estimation process

When the length of each bin is equal (which is very much often in the
literature), this approach has a single parameter, the number of bins, and so
it is possible to obtain a more refined or coarser description of the data by
controlling the number of bins. Therefore, it could be regarded as a
deterministic transformation from a NHPP to a HPP (for example, see 28,
Section 2.10) when the number of bins gets increased. In the learning process,
this approach usually restricts its attention to special parametric families
of rate functions to improve the learning quality. Learning on small
subintervals can make an accurate approximation for all observed data. The
powerful reason motivating the use of this model selection is the reduction of
the complexity of the chosen functions and the reduction of training errors.
Therefore, data binning is necessary and pervasive in both estimation and
learning process to reduce the empirical risk value. This could be shown by
the following proposition.
Given $m$ observations, the functionals in (3), (4) and (6) can be rewritten
as $R(\alpha,m)$, $R_{emp}(\alpha,m)$ and $\xi(m)$ related to the number of
observations, $m$. Suppose the range of considered time $[0,T]$ is divided
into $n$ consecutive nonoverlapping subintervals by knots
$b_{0}=0<b_{1}<\dots<b_{n}=T$. Let $t_{i}$, $0\leq t_{i}\leq T$, be arrival-
time point of event $x_{i}$, $S_{k}=\\{z_{i}=(x_{i},y_{i})\emph{ }|\emph{
}b_{k-1}<t_{i}\leq b_{k}\wedge i=\overline{1,m}\\}$, $\forall
k=\overline{1,n}$. We denote the empirical risk functional on the $k^{th}$
subinterval as $R^{k}_{emp}(\alpha_{k},m_{k})$, where, $m_{k}=|S_{k}|$. The
empirical risk functional after binning is computed as:

$\displaystyle\hat{R}_{emp}(\alpha_{1},\dots,\alpha_{n},m)=\frac{1}{m}\sum_{k=1}^{n}m_{k}R^{k}_{emp}(\alpha_{k},m_{k})$
(7)

######  Proposition 1.

For any $n\geq 2$, the inequality

$\displaystyle\min_{\alpha_{1},\dots,\alpha_{n}\in\Lambda}\hat{R}_{emp}(\alpha_{1},\dots,\alpha_{n},m)\leq\min_{\alpha\in\Lambda}R_{emp}(\alpha,m)$

holds true.

###### Proof.

For each subinterval $k$, let us denote
$\hat{\alpha}_{k}=\arg\min_{\alpha_{k}\in\Lambda}R^{k}_{emp}(\alpha_{k},m_{k})$
and $\hat{\alpha}=\arg\min_{\alpha\in\Lambda}R_{emp}(\alpha,m)$. The
inequality in the proposition could be inferred from the following inequality:

$\displaystyle R^{k}_{emp}(\hat{\alpha}_{k},m_{k})\leq
R_{emp}(\hat{\alpha},m_{k})$

for each subinterval $k=\overline{1,n}$. ∎

Proposition 1 shows that dividing the training data into smaller bins likely
leads to the decrease of empirical risk (better approximation of training
data).

In prior works, the common approach is to utilize the data binning with equal
length [6, 4, 29] and this is independent of the estimation step. The number
of subintervals is typically dependent on the modeler’s experience and not
data driven. Having a small number of subintervals will lead to a large
empirical risk value, while having a larger number of subintervals may improve
the empirical risk but may be risky. Here, one of the risks is that training
data is sparse and nonhomogeneous in time, therefore, there are almost always
empty subintervals. The other risk is generalization and the following
questions could be asked:

  * •

_Does the binning on finite data increase the probability of overfitting?_

  * •

_How to divide the data into bins to avoid overfitting?_

To answer the first question, according to [26], we consider the relationship
between a binning method on finite data and the overfitting problem. In this
context, we formulate a proposition to show that the generalization of a
learning machine depends not only on the empirical risk value and the VC
dimension of $F$ (see Section 3.4.4) but also on the binning method. In
detail, if we use the smaller number of observations (i.e., the smaller size
of subintervals) for approximation, the error on the training set can be
reduced but the error on the testing set may increase.
We rewrite the right-hand side of (5) in the new form

$\displaystyle\overline{R}(\alpha,m)=R_{emp}(\alpha,m)+2B\xi(m)\left(1+\sqrt{1+\frac{R_{emp}(\alpha,m)}{B\xi(m)}}\right)$

Then functional (5) can be written in as follow:

$\displaystyle R(\alpha,m)\leq\overline{R}(\alpha,m)$

######  Proposition 2.

There exists a number $m^{*}$ such that for all $m<m^{*}$, we have

$\displaystyle\overline{R}(\alpha,m)\geq\overline{R}(\alpha,m^{*})$

###### Proof.

Clearly, from the definition of the empirical risk functional, we have
$R_{emp}(\alpha,m)\xrightarrow[m\rightarrow\infty]{a.s}R(\alpha,m)$.
Therefore, the functional $R_{emp}(\alpha,m)$ is bounded from below. We denote

$\displaystyle m^{*}=\arg\min_{m}R_{emp}(\alpha,m),$

where $m>>h$ (see 26). Moreover, we have:

$\displaystyle\xi(m)=\frac{h[\ln{(2m/h)}+1]-\ln{(\eta/4)}}{m}$

where $0<\eta<1$. Thus, we obtain $\ln{\eta/4}<0$ and the functional
$(\ln{2x})/x$ is monotonically decreasing on $[e^{2}/2,+\infty)$. Therefore,
we obtain

$\displaystyle\xi(m)$
$\displaystyle\left(1+\sqrt{1+\frac{R_{emp}(\alpha,m)}{B\xi(m)}}\right)\geq$
$\displaystyle\xi(m^{*})\left(1+\sqrt{1+\frac{R_{emp}(\alpha,m^{*})}{B\xi(m^{*})}}\right)$

This proposition is proved. ∎

We see that the number of observations on each subinterval depends on the
number of subintervals. Basically, the larger number of subintervals leads to
a smaller number of observations on subintervals. Therefore, this proposition
specifies that binning the subintervals into smaller ones carries a risk of
overfitting.
To illustrate the quality of the learning process depending on the number of
subintervals for a tested problem and data described in section 6, Fig. 3
shows how the training (solid line) and test (dashed line) errors evolve when
the number of subintervals gets increased and the binning method is equal
length binning. We could see that the test error increases while the training
error steadily decreases starting from 15 subintervals (bins). A typical error
graph for the overfitting phenomenon is observed.

Figure 3: Overfitting in piecewise-polynomial regression on the San Francisco
taxi request data (see Section 6)

To answer the second question, we propose, in the next section, two
regularization methods for binning the training data, one is based on Ivanov
regularization and the other is based on Tikhonov regularization.

##  5 Learning Regularization of NHPPs

In this section, we propose a general framework for learning NHPPs. Two novel
binning methods inspired by Ivanov and Tikhonov regularizations are given.
Finally, the general regularized learning algorithm of NHPPs is detailed.

###  5.1 Learning Framework

In this part, we propose a process of learning the rate function of a NHPP.
Our process is decomposed into two steps: data binning and regression. In the
data binning step, samples are placed in discrete bins by a binning method.
After that, a regression learning technique is applied. Our learning process
is different from the literature in that the binning step is not conducted one
time but repeatedly with the regression step. This makes our binning step more
adaptive to the training data. The whole learning process is illustrated in
Fig. 4.

Figure 4: The proposed learning framework

###  5.2 Data Binning Methods

As previously discussed, binning is important for learning/ estimation of
NHPPs. While having too few bins could result in bad estimations of NHPPs
(underfitting), using too many bins could risk overfitting. In general, a
learning/estimation of NHPPs is depend on not only the number of bins but also
the way of dividing the data into bins. However, to the best of our knowledge,
all of the previous works on estimation/learning of NHPPs have sidelined this
problem and taken ad-hoc, non-data-driven approaches in choosing the way to
divide data into bins (e.g., divide data into equal-length intervals with a
predefined bin number). In this subsection, we propose two adaptive and data
driven binning methods that try to take into account the learning
regularization.

####  5.2.1 Ivanov regularization based data binning method

Based on the idea of Ivanov regularization, we propose a general data division
algorithm using distribution tests for breaking up the considered time
interval into subintervals. The distribution test plays the role of the
constraint in Ivanov regularization, which is used to restrict the bins to the
set of subintervals that respect the Poisson property of training data. The
details of the algorithm are shown in Algorithm 1. In the algorithm,
$dP_{best}$ indicates the best solution found during the search. It divides an
interval into subintervals (lines 2-3) by a uniformly random point (line 1).
Then, the Poisson property on each subinterval is tested using the Log Test
[6] (line 4). Lines 5-13 update the best training error and recursively work
on the newly created subintervals. Lines 14-18 continue the division if the
stopping criterion (for instance, the recursive depth or the total number of
bins) is not met.

Algorithm 1 Recursive data binning method based on the Ivanov regularization

1:dividedInterval$(l,u,O,dP)$

2:

3:l: lower bound of interval

4:u: upper bound of interval

5:O = $\\{t_{1},...,t_{m}\\}$: arrival-time points

6:dP: knots

7:minR, $dP_{best}$: the best-found risk value and knot list

8:$p\leftarrow$ random number in U $\simeq$ $\mathrm{Uniform}(l,u)$

9:$O_{1}\leftarrow\\{t_{i}\emph{ }|\emph{ }t_{i}\in O\wedge l\leq t_{i}<p\\}$

10:$O_{2}\leftarrow\\{t_{i}\emph{ }|\emph{ }t_{i}\in O\wedge p\leq t_{i}<u\\}$

11:Using a distribution test for testing Poisson property on $O_{1}$ and
$O_{2}$

12:if Both $O_{1}$ and $O_{2}$ are passed the distribution test then

13: $dP\leftarrow dP\cup{p}$

14: Let $er$ be the value of risk functional (7)

15: if $er\leq minR$ then

16: $minR\leftarrow er$

17: $dP_{best}\leftarrow dP$

18: end if

19: dividedInterval($l,p,O_{1},dP$)

20: dividedInterval($p,u,O_{2},dP$)

21:else

22: if stop-criterion not met then

23: dividedInterval($l,u,O,dP$)

24: end if

25:end if

26:

####  5.2.2 Data binning method based on the idea of Tikhonov regularization

Based on the idea of Tikhonov regularization, we propose a binning method
where an additional regularized term is added to the empirical risk function.
This term is used to control the number of subintervals based on their
lengths. The empirical risk functional (7) is reformulated as:

$\displaystyle\hat{R}_{emp}(\alpha_{1},\dots,\alpha_{n},m)=$
$\displaystyle\frac{1}{m}\Sigma_{k=1}^{n}m_{k}R^{k}_{emp}(\alpha_{k},m_{k})$
$\displaystyle+\gamma\Sigma_{k=1}^{n-1}\delta_{k}$ (8)

where $\delta_{k}=m_{k}R^{k}_{emp}(\alpha_{k},m_{k})/(b_{k}-b_{k-1})$ is the
risk value scaled by the length of subinterval $k^{th}$ found at the previous
division step. When we increase the number of subintervals, in general, the
length of subintervals and the risk value are decreased. However, the
decreasing rate of the risk value is usually slower than that of the
subinterval length. Thus, when the number of subintervals is large enough, the
penalty term tends to increase. This form of regularization term represents
the momentum of the empirical risk function related to the time interval after
each division iteration. The proposed idea is depicted in Algorithm 2. Line 1
generates a divided point randomly. Line 3 computes the new empirical risk
value. Lines 4-7 save the best solution found. Lines 8-11 continue the
division if the stopping criterion (for instance, the recursive depth or the
total number of bins) is not met.

Algorithm 2 Recursive data binning method based on the Tikhonov regularization

1:dividedInterval$(l,u,dP)$

2:

3:l: lower bound of interval

4:u: upper bound of interval

5:dP: knots

6:minR, $dP_{best}$: the best-found risk value and knot list

7:$p\leftarrow$ random number in U $\simeq$ $\mathrm{Uniform}(l,u)$

8:$dP\leftarrow dP\cup{p}$

9:Let $er$ be the value of risk functional (5.2.2)

10:if $er\leq minR$ then

11: $minR\leftarrow er$

12: $dP_{best}\leftarrow dP$

13:end if

14:if stop-criterion not met then

15: dividedInterval($l,p,dP$)

16: dividedInterval($p,u,dP)$)

17:end if

18:

###  5.3 A Proposed Algorithm for Estimating the Rate Function of NHPPs

Algorithm 3 is an iterative algorithm for estimating the rate function of a
NHPP. It repeatedly calls Algorithm 1 (or Algorithm 2) until the stopping
criterion is met. In the algorithm, $minR$ and $dP_{best}$ indicate the best
risk value and binning solution found during the search.

Algorithm 3 Algorithm for learning the NHPP rate function

1:dividedInterval

2:

3:[0, T]: considered time interval

4:O = $\\{t_{1},...,t_{m}\\}$: arrival-time points

5:minR, $dP_{best}$: the best-found risk value and knot list

6:Initialize $dP_{best}\leftarrow dP\leftarrow\\{\emptyset\\}$

7:$minR\leftarrow\infty$

8:while stop-criterion not met do

9: $dP\leftarrow\textsc{dividedInterval}(0,T,O,dP)$

10: Let $er$ be the value of risk functional (7)

11: if $er<minR$ then

12: $minR\leftarrow er$

13: $dP_{best}\leftarrow dP$

14: end if

15: $dP\leftarrow\\{\emptyset\\}$

16:end while

17:

##  6 Experiments and Results

###  6.1 Selection of Hypothesis Space

To test the effect of binning on learning generalization of NHPPs, we select
hypothesis spaces with fixed VC-dimensions in each experiment. This practice
is common in the literature of NHPP estimation [31, 9, 6, 7, 4]. The two most
popular function spaces for NHPP estimation in the literature are selected,
namely, the spaces of linear and polynomial functions [31, 4].

###  6.2 Tested Problems and Data

In the experiments, we have chosen some well-known NHPPs learning/estimation
problems from the literature that model different kinds of events such as
daily taxi call arrivals (e.g., 5, 32, 33) and telephone call arrivals at
customer service call centers (e.g., 6, 34, 35). For these problems, the true
rate functions are unknown and must be learned/ approximated from data.
Besides, we also consider the problem of learning a NHPP from a synthetic
dataset, which was randomly generated with a known piecewise-linear rate
function. These problems and their data are described as follows.

  * •

Problem 1 (AF): Learning the arrival rate of randomly generated requests
coming from a NHPP. The synthetic dataset is generated by a standard algorithm
described in [36] with the piecewise-linear rate function $\lambda(t)$ as
follows:

$\displaystyle\lambda(t)=\begin{cases}\cfrac{13}{36}\times t+7,&\quad\text{if
}0\leq t<10800,\\\ \cfrac{5}{36}\times t+15,&\quad\text{if }10800\leq
t<21600,\\\ \cfrac{25}{36}\times t-25,&\quad\text{if }21600\leq t<32400,\\\
-\cfrac{1}{2}\times t+104,&\quad\text{if }32400\leq t<43200,\\\
-\cfrac{1}{18}\times t+40,&\quad\text{if }43200\leq t<54000,\\\
\cfrac{1}{3}\times t-30,&\quad\text{if }54000\leq t<64800,\\\
-\cfrac{4}{9}\times t+138,&\quad\text{if }64800\leq t<75600,\\\
-\cfrac{5}{9}\times t+166,&\quad\text{if }75600\leq t<86400,\\\
0,\quad\text{otherwise.}\end{cases}$

Subintervals of lengths are equal, using known periodic arrival rate
functions. For each day, the number of requests is randomly chosen in [7000,
8000].

    * –

Training dataset: 365 days.

    * –

Testing dataset: 31 days.

  * •

Problem 2 (SF): Learning the arrival rate of San Francisco taxi requests.

    * –

Data: Dataset of mobility traces of taxi cabs in San Francisco, US. It
contains GPS coordinates of approximately 564 taxis collected over 31 days.
The average number of requests per day is around 14.255. The time calls are
specified to be the pickup time subtracting 10 minutes and clustered every 5
minutes.

    * –

Link for the dataset: _http://cabspotting.org/_

    * –

Training dataset: The data collected over 24 days in 03/2010.

    * –

Testing dataset: The data collected over 7 days in 03/2010.

  * •

Problem 3 (PT): Learning the arrival rate of Porto taxi requests.

    * –

Data: An accurate dataset describes a complete year of the trajectories for
all of the 442 taxis running in the city of Porto, Portugal. These taxis
operate through a taxi dispatch central and serve 1710589 trips in a year from
01/07/2013.

    * –

Link for the dataset: _https://www.kaggle.com_

    * –

Training dataset: The data collected from from 01/07/2013 to 30/05/2014.

    * –

Testing dataset: The data collected over 30 days in 06/2014.

  * •

Problem 4 (IB): Learning the arrival rate of calls to a call-center at a small
Israeli bank.

    * –

Data: detailed description of the data could be found in [37]. There are
135780 calls arrived in four months. Each record includes the time stamp for
when the call arrived. The average number of arrival calls per minute is 1.1.

    * –

Link for the dataset: _http://ie.technion.ac.il/serveng/callcenterdata/inde_
-_x.html_

    * –

Training dataset: The data collected from from 01/01/1999 to 30/04/1999.

    * –

Testing dataset: The data collected over 30 days in 05/1999.

In this paper, we conduct three main experiments.

####  6.2.1 Experiment 1

  * •

_Description_ : The first experiment is to confirm that the learning with data
binning on finite data can carry the probability of being overfitting.

  * •

_Data_ : The SF, PT, IB, and AF datasets are used.

####  6.2.2 Experiment 2

  * •

_Description_ : In the second experiment, we perform our proposed algorithms
to evaluate the efficiency of learning regularization for NHPPs.

  * •

_Data_ : The SF, PT, IB and AF datasets are used.

####  6.2.3 Experiment 3

  * •

_Description_ : In the last experiment, we test our proposed algorithms for a
more general case in real life, i.e., Dynamic Vehicle Routing Problem (DVRP)
with spatio-temporal data.

  * •

_Data_ : The SF dataset is used.

###  6.3 Training and Testing Processes

In our experiments, the range of considered time $[0,T]$ is set by day unit.
The arrival rate functions depend on the time of day (i.e., we assume that
requests of different dates come from the same distribution). Therefore, the
number of arrival events is a result of repeated day-to-day observations. For
each day in the training dataset, the number of requests coming at a time of
day (the arrival time is rounded by minute) is considered as a data point. To
avoid an ecological fallacy, all those data points are used to train the model
(see 38). As illustrated in Fig. 5, the number of values on the vertical axis
(possibly equal) at each value on the horizontal axis is equal to the number
of days in the considered dataset. The training process minimizes the _Root
Mean Squared Error_ (RMSE) on all those data points. The testing process is
performed similarly on the test dataset.

Figure 5: An example of the training and testing processes

###  6.4 Experiment 1: Overfitting Problem

To show the possibility of being overfitting when using data binning, we
conducted the experiments of learning NHPPs with increasing number of bins. It
could simply be done by relaxing the constraint that requires the Poisson
property in Algorithm 1. In line 6, if either of the two datasets
$O_{1},O_{2}$ does not pass the distribution test, and its length is not
bigger than $2\eta$ (minute unit), the division on dataset $O$ is not
performed. This ensures that the number of bins is finite as common in the
literature. The piecewise-polynomial regressions with different $\eta$ were
performed. The RMSE training ($\sigma_{1}$) and testing ($\sigma_{2}$) errors
for all tested problems are presented in Fig. 6. Table 1 shows the detailed
numerical results on the SF dataset.

Table 1: Results of the piecewise-polynomial regression with the relaxed division algorithm. Ins |  $\\#$b | $\sigma_{1}$ | $\sigma_{2}$
---|---|---|---
$\eta$ = 600 | 3 | 6.54 | 6.3
$\eta$ = 480 | 4 | 6.43 | 6.23
$\eta$ = 120 | 7 | 5.23 | 5.23
$\eta$ = 100 | 8 | 5.21 | 5.2
$\eta$ = 80 | 11 | 5.1 | 5.07
$\eta$ = 60 | 14 | 4.91 | 5.21
$\eta$ = 50 | 15 | 4.9 | 5.19
$\eta$ = 40 | 19 | 4.88 | 5.21
$\eta$ = 30 | 26 | 4.76 | 5.24
$\eta$ = 20 | 48 | 4.38 | 5.31
$\eta$ = 10 | 79 | 4.26 | 5.48
Figure 6: Overfitting problem in piecewise-polynomial regression

It can be seen that when the number of subintervals ($\\#$b) gets increased,
initially, both $\sigma_{1}$ and $\sigma_{2}$ errors steadily decrease.
However, starting from 14 bins, the training error decreases slowly and the
test error gets increased faster. This behavior is commonly attributed to the
overfitting process in the literature of machine learning. The same is true
for all other tested problems (see Table 5 in the Appendix) and it is
consistent with the theoretical result shown in subsection 4.2.

###  6.5 Experiment 2: Learning Regularization of NHPPs

We denote the RMSE on the training and the test datasets as $\sigma_{1}$ and
$\sigma_{2}$, respectively. The number of subintervals found by Algorithm 3 is
denoted as $\\#b$. For comparison, both regularized data binning method based
on Ivanov ($\textit{DBM}_{I}$) and Tikhonov ($\textit{DBM}_{T}$)
regularization were performed. We also conducted the common data binning
method in the literature, where the considered interval is divided into $\\#b$
subintervals with equal length. $\overline{\sigma_{1}}$ and
$\overline{\sigma_{2}}$ denote the RMSEs on the training and test datasets
found by the estimation process. Notation $\rho$ presents the improvement
percentage of $\sigma_{2}$ found by the regression on subintervals with our
binning method compared to $\overline{\sigma_{2}}$ found by the regression on
subintervals with equal length.

####  6.5.1 The piecewise-linear and linear regression

Table 2 shows the results of learning by linear (_Unbinned_ column) and
piecewise-linear functions. It is undoubted that binning is essentially
important for learning NHPPs with both training and test errors are
significantly lower on all tested problems as shown in the table. The results
also show that adaptive binning with regularization helps to improve on test
errors from 10.73% to 25.57% compared to equal length binning. For the two
regularization methods, Ivanov regularization ($\textit{DBM}_{I}$) gives
better results than Tikhonov regularization ($\textit{DBM}_{T}$) on the AF, PT
and IB problems and worse on SF. The result illustrations are placed in the
Appendix for the convenience of reference (Figures 8-15).

Table 2: The regression results of the rate function approximated by linear and piecewise-linear functions. Ins | Unbinned | Binning with $\textit{DBM}_{I}$ | Binning with $\textit{DBM}_{T}$
---|---|---|---
$\sigma_{1}$ | $\sigma_{2}$ | $\sigma_{1}$ | $\sigma_{2}$ | $\\#b$ | $\overline{\sigma_{1}}$ | $\overline{\sigma_{2}}$ | $\rho$ | $\sigma_{1}$ | $\sigma_{2}$ | $\\#b$ | $\overline{\sigma_{1}}$ | $\overline{\sigma_{2}}$ | $\rho$
AF | 8.89 | 10.6 | 4.18 | 4.74 | 10 | 5.17 | 5.31 | 10.73 | 4.75 | 4.86 | 12 | 5.21 | 5.46 | 10.98
SF | 11.1 | 13.4 | 3.2 | 4.27 | 12 | 5.02 | 5.14 | 16.93 | 3.16 | 3.87 | 13 | 4.96 | 5.2 | 25.57
PT | 5.43 | 5.69 | 2.82 | 3.28 | 15 | 4.73 | 4.05 | 19.01 | 2.96 | 3.84 | 11 | 4.92 | 4.47 | 14.09
IB | 5.74 | 6.06 | 2.42 | 2.73 | 11 | 2.78 | 3.19 | 14.42 | 2.41 | 2.81 | 9 | 3.12 | 3.36 | 16.36

####  6.5.2 The polynomial and piecewise-polynomial regression

In this part, the hypothesis spaces are polynomial (unbinned) and piecewise-
polynomial functions (binning). Table 3 presents the quality improvement of
algorithms using regularizations $\textit{DBM}_{I}$ and $\textit{DBM}_{T}$.

Table 3: The regression results of the rate function approximated by polynomial and piecewise-polynomial functions. Ins | Unbinned | Binning with $\textit{DBM}_{I}$ | Binning with $\textit{DBM}_{T}$
---|---|---|---
$\sigma_{1}$ | $\sigma_{2}$ | $\sigma_{1}$ | $\sigma_{2}$ |  $\\#$b | $\overline{\sigma_{1}}$ | $\overline{\sigma_{2}}$ | $\rho$ | $\sigma_{1}$ | $\sigma_{2}$ |  $\\#$b | $\overline{\sigma_{1}}$ | $\overline{\sigma_{2}}$ | $\rho$
AF | 5.8 | 5.97 | 4.84 | 4.91 | 8 | 5.33 | 5.38 | 8.74 | 4.78 | 4.96 | 12 | 5.27 | 5.49 | 9.65
SF | 3.69 | 4.41 | 3.02 | 3.56 | 15 | 3.34 | 4.33 | 17.78 | 2.34 | 3.69 | 14 | 3.41 | 4.38 | 15.75
PT | 3.15 | 4.72 | 1.91 | 3.9 | 20 | 3.02 | 4.62 | 15.58 | 1.91 | 3.98 | 12 | 2.88 | 4.53 | 12.14
IB | 2.39 | 3.38 | 1.98 | 2.63 | 16 | 2.17 | 3.12 | 15.71 | 2.06 | 2.69 | 13 | 2.24 | 2.84 | 5.28

The results again confirm the importance of data binning and the
regularization approaches. $\textit{DBM}_{I}$ has the best test errors on all
tested problems improving over those of equal length binning method from 8.74%
to 17.78%. Comparing with the results in Table 2, the piecewise-polynomial
regression on PT and AF datasets are worse the piecewise-linear regression. A
natural explanation is that the rate functions of these datasets come from
piecewise-linear functions (certainly for AF) or the choice of the degree of
polynomial (PT).

Overall, the experimental results clearly indicate the effectiveness of
regularization methods over the vanilla equal length binning method. Of the
two methods, $\textit{DBM}_{I}$ seems to be the best. However,
$\textit{DBM}_{T}$ might still be useful since it has lower computational time
(without repeatedly performing statistical tests) and for the case of high
dimensional NHPPs, where the efficient statistical tests are not available.

###  6.6 Experiment 3: Dynamic Vehicle Routing Problem (DVRP) with Spatio-
Temporal Data

Figure 7: The result on the SF dataset with 20 sub-areas Table 4: The profit of scheduling algorithm using our proposed learning method. Ins | $p_{1}$ | $p_{2}$ | $\rho(\%)$
---|---|---|---
Day 1 | 27.44 | 29.14 | 6.19
Day 2 | 30.39 | 32.11 | 5.66
Day 3 | 31.71 | 33.98 | 7.16
Day 4 | 29.62 | 31.03 | 4.76
Day 5 | 29.8 | 31.61 | 6.07
Day 6 | 29.24 | 30.15 | 3.11
Day 7 | 29.62 | 31.17 | 5.23
Day 8 | 28.56 | 29.72 | 4.06
Day 9 | 10.16 | 10.85 | 6.79

In this experiment, we evaluate the benefit of using our regularized learning
of NHPPs for a recent DVRP application. In the DVRP, a spatio-temporal process
can be reduced to temporal processes by a clustering algorithm as proposed in
[39]. A large geographical area can be clustered into non overlap sub-areas by
the $K$-means clustering algorithm. Most of recent works on DVRP have utilized
this approach, for instance, [29, 5]. In [5], the authors consider the share-
a-ride taxi scheduling problem and assume that the arrival rate of taxi
requests only depend on period and area and hence are fixed overall days. The
city is partitioned into a set of sub-areas. Each day is binned into
subintervals with equal lengths. The number of arrival requests on each sub-
area in a subinterval is generated by a Poisson distribution. In this
experiment, we change the learning method in [5] by our proposed learning
process using the regularization method $DBM_{I}$ to learn the arrival rate of
taxi requests. In our experiment, the arrival rate of taxi requests is learned
and used to suggest the best route from the last drop-off point to a parking
lot for each taxi. Fig. 7 relates the piecewise-polynomial regression result
on the SF dataset in which the number of sub-areas is equal to 20 ($K=20$).
Table 4 shows the profit of scheduling algorithms. We denote the OTSF-DP
algorithm in [5] by $p_{1}$ and the OTSF-DP using our new predicted rate
function by $p_{2}$. The percentage of improvement is denoted by
$\rho=(p_{2}-p_{1})100/p_{1}$. It can be seen that by using the new predicted
rate function obtained via regularized learning approach in the paper, the
profit is improved significantly for all nine tested days. The percentage of
improvement is from 3.11% (909 USD per day) to 7.16% (2270 USD per day), which
is of significantly economical values.

##  7 Conclusion

In this paper, we formulate the problem of estimation of NHPPs from finite and
limited data as a machine learning problem. We have theoretically and
experimentally shown that while the binning of observed data is essential for
learning/estimating NHPPs, it could lead to overfitting. Consequently, we have
proposed a learning framework of NHPPs that utilizes two regularization
methods based on ideas of the Tikhonov [27] and Ivanov regularization [40]. We
tested our regularization learning methods on synthetic and real-world
datasets and the experimental results show the effectiveness of these methods.
The learning methods were further tested by embedding their results in a
recent real-world application of share-a-ride taxi scheduling using spatio-
temporal data. Our new regularized learning methods help to significantly
improve the performance of one of the state-of-the-art dynamic scheduling
algorithms producing better economical profits. Moreover, the new regularized
learning methods completely remove the need of ad-hoc parameter tuning in
binning (such as the number of bins) by a data driven binning method.

In future, we will extend our approach to learning high dimensional NHPPs and
extend the VC-theory to fully cover the statistical learning of NHPPs.

## Declarations

  * •

Funding: The authors did not receive support from any organization for the
submitted work

  * •

Conflict of interest/Competing interests: The authors have no competing
interests to declare that are relevant to the content of this article.

  * •

Ethics approval: Not applicable

  * •

Consent to participate: Not applicable

  * •

Consent for publication: Not applicable

  * •

Availability of data and materials: The instances and our results are readily
available on _https://github.com/sonnv188/poissonprocess.git_

  * •

Code availability: The instances and our results are readily available on
_https://github.com/sonnv188/poissonprocess.git_

  * •

Authors’ contributions: All authors contributed to the study conception and
design. Conceptualization, methodology and validation was proposed by Hoai
Nguyen Xuan. Material preparation, implementing and data analysis were
performed by Son Nguyen Van. The first draft of the manuscript was written by
Son Nguyen Van and all authors commented on previous versions of the
manuscript. All authors read and approved the final manuscript.

##  Appendix A Learning Results

Table 5: Results of the piecewise-polynomial regression with the relaxed division algorithm Ins | AF | PT | IB
---|---|---|---
2-10 |  $\\#$b | $\sigma_{1}$ | $\sigma_{2}$ |  $\\#$b | $\sigma_{1}$ | $\sigma_{2}$ |  $\\#$b | $\sigma_{1}$ | $\sigma_{2}$
$\eta$ = 600 | 2 | 6.7 | 5.23 | 2 | 5.84 | 5.18 | 3 | 3.78 | 3.28
$\eta$ = 480 | 4 | 6.57 | 5.16 | 3 | 5.81 | 5.05 | 6 | 3.65 | 3.15
$\eta$ = 120 | 7 | 6.59 | 4.95 | 7 | 5.57 | 4.92 | 8 | 3.49 | 3.01
$\eta$ = 100 | 9 | 6.28 | 4.9 | 9 | 5.32 | 4.76 | 10 | 3.25 | 2.89
$\eta$ = 80 | 11 | 5.91 | 5.08 | 11 | 5.26 | 4.53 | 11 | 3.16 | 2.74
$\eta$ = 60 | 13 | 5.75 | 5.11 | 14 | 5.03 | 4.17 | 14 | 3.02 | 2.9
$\eta$ = 50 | 17 | 5.6 | 5.18 | 17 | 4.98 | 3.92 | 17 | 2.87 | 3.13
$\eta$ = 40 | 21 | 5.41 | 5.18 | 19 | 4.85 | 3.95 | 20 | 2.73 | 3.25
$\eta$ = 30 | 26 | 5.34 | 5.46 | 28 | 4.8 | 4.01 | 27 | 2.66 | 3.25
$\eta$ = 20 | 44 | 5.22 | 5.61 | 39 | 4.21 | 4.57 | 41 | 2.62 | 3.32
$\eta$ = 10 | 85 | 5.09 | 5.82 | 82 | 4.14 | 4.59 | 83 | 2.54 | 3.39

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 8: The piecewise-linear regression of the arrival rate function on the
AF dataset

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 9: The piecewise-linear regression of the arrival rate function on the
SF dataset

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 10: The piecewise-linear regression of the arrival rate function on the
PT dataset

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 11: The piecewise-linear regression of the arrival rate function on the
IB dataset

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 12: The piecewise-polynomial regression of the arrival rate function on
the AF dataset

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 13: The piecewise-polynomial regression of the arrival rate function on
the SF dataset

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 14: The piecewise-polynomial regression of the arrival rate function on
the PT dataset

(a) $\textit{DBM}_{I}$

(b) $\textit{DBM}_{T}$

Figure 15: The piecewise-polynomial regression of the arrival rate function on
the IB dataset

## References

  *   * Komaki [2021] Komaki, F.: Shrinkage priors for nonparametric bayesian prediction of nonhomogeneous poisson processes. IEEE Transactions on Information Theory 67(8), 5305–5317 (2021) https://doi.org/10.1109/TIT.2021.3084062
  * Jahani et al. [2021] Jahani, S., Zhou, S., Veeramani, D., Schmidt, J.: Multioutput gaussian process modulated poisson processes for event prediction. IEEE Transactions on Reliability 70(4), 1569–1580 (2021) https://doi.org/10.1109/TR.2021.3088094
  * Moore and Davenport [2019] Moore, M.G., Davenport, M.A.: Estimation of poisson arrival processes under linear models. IEEE Transactions on Information Theory 65(6), 3555–3564 (2019) https://doi.org/10.1109/TIT.2018.2889489
  * Alizadeh and Papp [2013] Alizadeh, F., Papp, D.: Estimating arrival rate of nonhomogeneous poisson processes with semidefinite programming. Annals of Operations Research 208(1), 291–308 (2013) https://doi.org/10.1007/s10479-011-1020-2
  * Son et al. [2022] Son, N.V., Nhan, V.T.H., Dung, P.Q., Hoai, N.X., Babaki, B., Dries, A.: Novel online routing algorithms for smart people-parcel taxi sharing services. ETRI Journal 44, 208–219 (2022) https://doi.org/10.4218/etrij.2021-0406
  * Brown et al. [2005] Brown, L., Gans, N., Mandelbaum, A., Sakov, A., Shen, H., Zeltyn, S., Zhao, L.: Statistical analysis of a telephone call center. Journal of the American Statistical Association 100(469), 36–50 (2005) https://doi.org/10.1198/016214504000001808
  * Alizadeh et al. [2008] Alizadeh, F., Eckstein, J., Noyan, N., Rudolf, G.: Arrival rate approximation by nonnegative cubic splines. Operations Research 56(1), 140–156 (2008) https://doi.org/10.1287/opre.1070.0443
  * Avanzi et al. [2021] Avanzi, B., Taylor, G., Wong, B., Xian, A.: Modelling and understanding count processes through a markov-modulated non-homogeneous poisson process framework. European Journal of Operational Research 290(1), 177–195 (2021) https://doi.org/10.1016/j.ejor.2020.07.022
  * Leemis [2004] Leemis, L.M.: Technical note: Nonparametric estimation and variate generation for a nonhomogeneous poisson process from event count data. IIE Transactions 36(12), 1155–1160 (2004) https://doi.org/10.1080/07408170490507693
  * Arkin and Leemis [2000] Arkin, B.L., Leemis, L.M.: Nonparametric estimation of the cumulative intensity function for a nonhomogeneous poisson process from overlapping realizations. Management Science 46(7), 989–998 (2000) https://doi.org/10.1287/mnsc.46.7.989.12037
  * Leemis [1991] Leemis, L.M.: Nonparametric estimation of the cumulative intensity function for a nonhomogeneous poisson process. Management Science 37(7), 886–900 (1991) https://doi.org/10.1287/mnsc.37.7.886
  * Schatz et al. [2022] Schatz, M., Wheatley, S., Sornette, D.: The arma point process and its estimation. Econometrics and Statistics 24, 164–182 (2022) https://doi.org/10.1016/j.ecosta.2021.11.002
  * Kuhl and Bhairgond [2000] Kuhl, M.e., Bhairgond, P.s.: Nonparametric estimation of nonhomogeneous poisson processes using wavelets. 2000 Winter Simulation Conference Proceedings (Cat. No.00CH37165) (2000) https://doi.org/10.1109/wsc.2000.899764
  * Walter and Shen [1998] Walter, G.G., Shen, X.: Positive estimation with wavelets. Wavelets, Multiwavelets, and Their Applications Contemporary Mathematics, 63–79 (1998) https://doi.org/10.1090/conm/216/02964
  * Papp [2017] Papp, D.: Univariate polynomial optimization with sum-of-squares interpolants. Modeling and Optimization: Theory and Applications Springer Proceedings in Mathematics & Statistics, 143–162 (2017) https://doi.org/10.1007/978-3-319-66616-7_9
  * Poisson [1837] Poisson, S.D.: Probabilitédes jugements en matière criminelle et en matière civile, précédées des règles générales du calcul des probabilitiés. Bachelier (1837)
  * Suhov and Kelbert [2005] Suhov, Y.M., Kelbert, M.: Probability and statistics by example. Cambridge University Press (2005)
  * Applebaum [2004] Applebaum, D.: Lévy processes and stochastic calculus. Cambridge University Press (2004)
  * Ross [2016] Ross, S.M.: Stochastic processes. Wiley India (2016)
  * Sen and Cinlar [1975] Sen, P.K., Cinlar, E.: Introduction to stochastic processes. Journal of the American Statistical Association 70(352), 967 (1975) https://doi.org/10.2307/2285487
  * Lilliefors [1967] Lilliefors, H.W.: On the kolmogorov-smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association 62(318), 399–402 (1967) https://doi.org/10.1080/01621459.1967.10482916
  * Durbin [1973] Durbin, J.: Distribution theory for tests based on the sample distribution function. Society for Industrial and Applied Mathematics (1973) https://doi.org/10.1137/1.9781611970586
  * Sachs [1997] Sachs, L.: Angewandte statistik: Anwendung statistischer methoden. Springer (1997)
  * Kim and Whitt [2014] Kim, S.-H., Whitt, W.: Choosing arrival process models for service systems: Tests of a nonhomogeneous poisson process. Naval Research Logistics (NRL) 61(1), 66–90 (2014) https://doi.org/10.1002/nav.21568
  * Vapnik [1998] Vapnik, V.N.: Statistical learning theory. Wiley (1998)
  * Vapnik [2000] Vapnik, V.N.: Direct methods in statistical learning theory. The Nature of Statistical Learning Theory, 225–265 (2000) https://doi.org/%****␣Manuscript.bbl␣Line␣400␣****10.1007/978-1-4757-3264-1_8
  * Tikhonov [1963] Tikhonov, A.N.: On solving ill-posed problem and method of regularization. Dok/. Akad. Nauk USSR 153 70, 501–504 (1963)
  * Streit [2010] Streit, R.L.: Poisson point processes imaging, tracking, and sensing. Springer US (2010)
  * Ferrucci [2013] Ferrucci, F.: Pro-active dynamic vehicle routing: Real-time control and request-forecasting approaches to improve customer service. Physica-Verlag, 1–319 (2013)
  * Hsieh et al. [2017] Hsieh, K., Harlap, A., Vijaykumar, N., Konomis, D., Ganger, G., Gibbons, P.B., Mutlu, O.: Gaia: Geo-distributed machine learning approaching lan speeds. In: NSDI (2017)
  * Massey et al. [1996] Massey, W.A., Parker, G.A., Whitt, W.: Estimating the parameters of a nonhomogeneous poisson process with linear rate. Telecommunication Systems 5(2), 361–388 (1996) https://doi.org/10.1007/bf02112523
  * Yoon et al. [2019] Yoon, T., Cherry, C.R., Ryerson, M.S., Bell, J.E.: Carsharing demand estimation and fleet simulation with ev adoption. Journal of Cleaner Production 206, 1051–1058 (2019) https://doi.org/10.1016/j.jclepro.2018.09.124
  * Liu et al. [2019] Liu, L., Qiu, Z., Li, G., Wang, Q., Ouyang, W., Lin, L.: Contextualized spatial–temporal network for taxi origin-destination demand prediction. IEEE Transactions on Intelligent Transportation Systems 20(10), 3875–3887 (2019) https://doi.org/10.1109/TITS.2019.2915525
  * Shen and Huang [2008] Shen, H., Huang, J.Z.: Forecasting time series of inhomogeneous Poisson processes with application to call center workforce management. The Annals of Applied Statistics 2(2), 601–623 (2008) https://doi.org/10.1214/08-AOAS164
  * Xu et al. [2017] Xu, L., Kwan, M.-P., McLafferty, S., Wang, S.: Predicting demand for 311 non-emergency municipal services: An adaptive space-time kernel approach. Applied Geography 89, 133–141 (2017) https://doi.org/10.1016/j.apgeog.2017.10.012
  * Ross [2013] Ross, S.: Generating continuous random variables. Simulation, 69–96 (2013) https://doi.org/10.1016/b978-0-12-415825-2.00005-x
  * Mandelbaum et al. [2000] Mandelbaum, A., Sakov, A., Zeltyn, S.: Empirical analysis of a call center. Technical report, Technion, Israel Institute of Technology (January 2000)
  * Garrett [2002] Garrett, T.: Aggregated versus disaggregated data in regression analysis: Implications for inference. Economics Letters 81, 61–65 (2002) https://doi.org/10.1016/S0165-1765(03)00149-6
  * Bent and Hentenryck [2009] Bent, R.W., Hentenryck, P.V.: Dynamic vehicle routing with stochastic requests. In: International Joint Conference on Artificial Intelligence (2009). https://api.semanticscholar.org/CorpusID:60917372
  * Ivanov and Gupta [1978] Ivanov, V.V., Gupta, G.D.: The theory of approximate methods and their application to the numerical solution of singular integral equations. Journal of Applied Mechanics 45(1), 230–230 (1978) https://doi.org/10.1115/1.3424255