Title: INFERENCE OF UNOBSERVED EVENT STREAMS WITH NEURAL HAWKES PARTICLE SMOOTHING

Abstract: Events that we observe in the world may be caused by other, unobserved events. We consider sequences of discrete events in continuous time. When only some of the events are observed, we propose particle smoothing to infer the missing events. Particle smoothing is an extension of particle filtering in which proposed events are conditioned on the future as well as the past. For our setting, we develop a novel proposal distribution that is a type of continuous-time bidirectional LSTM. We use the sampled particles in an approximate minimum Bayes risk decoder that outputs a single low-risk prediction of the missing events. We experiment in multiple synthetic and real domains, modeling the complete sequences in each domain with a neural Hawkes process (Mei & Eisner, 2017) . On held-out incomplete sequences, our method is effective at inferring the ground-truth unobserved events. In particular, particle smoothing consistently improves upon particle filtering, showing the benefit of training a bidirectional proposal distribution.

1. INTRODUCTION:

Event streams, i.e., discrete events in continuous time, are often partially observed in the world. Given trained models of the complete data and the missingness mechanism, one can probabilistically predict the unobserved events, no matter whether they are missing at random (MAR) or missing not at random (MNAR) . Such an ability is useful in many applied domains:

• Scientific experiments. Scientific findings rely heavily on experimental observations, but due to practical limitations, measurements are often incomplete. For example, a systems neuroscientist will miss all the neural spikes from unmonitored neurons, but these are imputable to some degree from spikes observed elsewhere in the system, providing a more complete picture of neural activity. • Medical records. Some patients keep a diary or use a smartphone app to log their behavior in between their hospital visits, but other patients don't. Out-of-hospital events such as symptoms, self-administered medications, diet, and sleep may correlate with monitored inhospital events such as tests, diagnoses, and treatments. Inferring these events when they are not recorded may help doctors counsel patients. and "what I saw he did" would help the player make good decisions. Similar remarks apply to practical scenarios (e.g., military) where multiple actors compete and/or cooperate. • User interface interactions. Many important user actions are not observed by the software.

For example, users of an online news provider may have known what they need by reading an information-rich headline without clicking (to read more details). Such events are expensive to observe (e.g. via a gaze tracker) thus often missing. However, imputing them (to some extent) given the observed events (e.g. clicks) would be useful, e.g., helping the software to improve the user satisfaction in the long run. • Other partially observed event streams arise in online shopping, social media, etc.

A flexible probabilistic model for complete event streams is the neural Hawkes process (Mei & Eisner, 2017) , a recurrent neural model that allows past events to have complex influences on the x z 0 t (b) Particle smoothing Figure 1 : Stochastically imputing a taxi's pick-up events ( ) given its observed drop-off events ( ). At this stage of the computation, we are trying to determine the next event after the at time t1.

In Figure 1a (particle filtering), the neural Hawkes process's LSTM has already read the proposed and observed events at times ≤ t1. Its resulting state determines the model intensities and of the two event types (left image). Both intensities are low because there are few passengers at this time of day, so it happens that no event is proposed in (t1, t2). Specifically, Algorithm 1 determines that the next proposed event ( in middle image) would be somewhere after t2, though without bothering to determine its time precisely (line 36). Thus, is discarded (line 37), having been preempted by the observed @t2. We continue to extend the particle by feeding @t2 into the LSTM (right image) and proposing subsequent events based on the new intensities after t2. But because was low at t2, the @t2 was unexpected, which results in downweighting the particle (line 28): intuitively, we have just realized that this particle will be improbable under the posterior, because its complete sequence will include consecutive drop-offs far apart in time ( @t1, @t2). Figure 1b (particle smoothing) samples from a better-informed proposal distribution. A second LSTM (Appendix B) reads the future observations from right to left. Its state is used together with to determine the proposal intensities and (left image). Since a drop-off at t2 strongly suggests a pick-up before t2, considering the future increases the intensity of pick-up on (t1, t2) from to (while decreasing that of drop-off from to ). Consequently, the next proposed event is more likely to be a pick-up in (t1, t2) than it was in Figure 1a . If so (middle image), this event @t1,1 is fed into the neural Hawkes process's LSTM (right image). The updated state determines the new model intensities and : the higher intensity says that a dropoff is now expected, so we will not downweight the particle as much if t2 is the next event.

The updated also combines with to determine the new proposal intensities and , which are used to sample the next event (either an unobserved event at t1,2 ∈ (t1, t2), or the next observed event t2).

type and timing of subsequent events. Suppose we have a trained version of this model, as well as a separate probability model (the "missingness mechanism") that stochastically determines which of the events will be observed. We can then easily use Bayes' Theorem to define the posterior distribution over the complete sequence x z given just the observed events x: see equation (3) below. However, it is computationally difficult to reason about this posterior distribution. In this paper we provide approximate methods for sampling complete sequences from the posterior. In the remainder of the introduction, we provide a high-level sketch of these methods. Mei & Eisner (2017) give an algorithm to sample a complete sequence from a neural Hawkes process. Each event in turn is sampled given the complete history of previous events. However, this algorithm only samples from the prior. We will adapt it into a particle filtering algorithm that samples from the posterior, i.e., given all the observed events. The basic idea (Figure 1a ) is to draw the events in sequence as before, but now we force any observed events to be "drawn" at the appropriate times. That is, we add the observed events to the sequence as they happen (and they duly affect the distribution of subsequent events). There is an associated cost: if we are forced to draw an observed event that is improbable given its past history, we must downweight the resulting complete sequence accordingly, because evidently this particular past history-as determined by the previous drawswas inconsistent with the observed event, and hence cannot be part of a high-likelihood complete sequence. Using this method, we sample many sequences (or particles) of different relative weights.

Alas, this approach is computationally inefficient. Sampling a complete sequence that is actually probable under the posterior requires great luck, as the thinning algorithm must have the good fortune to draw only events that happen to be consistent with future observations. Such lucky particles would appropriately get a high weight relative to other particles. The problem is that we will rarely get such particles at all (unless we sample a very large number of particles). To get a more accurate picture of the posterior, we can improve the method by drawing each event from a smarter distribution that is explicitly conditioned on the future observations (rather than drawing the event in ignorance of the future and then downweighting the particle if the future does not turn out as hoped). This idea is called particle smoothing (Doucet & Johansen, 2009) . How does it work in our setting? The neural Hawkes process defines the distribution of the next event using the state of a continuoustime LSTM that has read the past history from left-to-right. When sampling a proposed event, we now use an modified distribution (Figure 1b ) that also considers the state of a second continuoustime LSTM that has read the future observations from right-to-left. As this augmented distribution is imperfect-merely a proposal distribution-we still have to reweight our particles to match the actual posterior under the model. But this reweighting is not as drastic as for particle filtering, because the new proposal distribution was constructed and trained to resemble the actual posterior.

Bidirectional recurrent neural networks have proven effective at predicting linguistic words and their properties given their left and right contexts (Graves et al., 2013; Bahdanau et al., 2015; Peters et al., 2018) : in particular, Lin & Eisner (2018) recently applied them to particle smoothing for sequence tagging, To our knowledge, however, this is the first time such an architecture has been extended to predict events in continuous time.

2. PRELIMINARIES 1:

We consider a missing-data setting (Little & Rubin, 1987) . We are given a fixed time interval [0, T ] over which events can be observed. Each possible outcome in our probability distributions is a complete event sequence in which each event is designated as either observed or missing. The random variables Obs, Miss, and Comp refer respectively to the sets of observed events, missing events, and all events over [0, T ]. Thus Comp = Obs Miss, where denotes disjoint union. Under the probability distributions we will consider, |Comp| is almost surely finite.

We will observe Obs to be some particular set of events x = {k 0 @t 0 , k 1 @t 1 , k 2 @t 2 , . . . , k I @t I , k I+1 @t I+1 }, where each k i ∈ {1, 2, . . . , K} is an event type and 0 = t 0 < t 1 < t 2 < . . . < t I < t I+1 = T are the times of occurrence in increasing order. This sequence always includes the special boundary events k 0 @t 0 = BOS@0 ("beginning of sequence") and k I+1 @t I+1 = EOS@T ("end of sequence").

Following each observed event k i @t i for 0 ≤ i ≤ I, we may hypothesize some J i ≥ 0 unobserved events {k i,1 @t i,1 , k i,2 @t i,2 , . . . , k i,Ji @t i,Ji }, where t i < t i,1 < t i,2 < . . . < t i,Ji < t i+1 . We may regard the subscript i as a shorthand for (i, 0), so the observed events are x = {k i,j @t i,j : j = 0} and the hypothesized unobserved (missing) events are z = {k i,j @t i,j : j = 0}.

The hypothesized complete event stream x z is thus indexed by pairs = i, j ordered lexicographically, where t < t if < . We write + 1 for the index of the event immediately after , so i, j + 1 def = i, j + 1 when j < J i and i, j + 1 def = i + 1, 0 otherwise. Each event k @t (whether observed or not) may be influenced by the history H(t )-the set of all observed and unobserved events that precede it, where we define H(t) def = {k @t : t < t} for any t ∈ [0, T ].

We assume that the complete sequence Comp was first generated by a complete data model p, and then some of the events were censored (marked as missing) by a missingness mechanism p miss :

p(Obs = x, Miss = z) = p(Comp = x z) • p miss (Miss = z | Comp = x z)

(1) Furthermore, our complete data model (such as a neural Hawkes process) will take the factored form

EQUATION

where p(k @t | H(t )) is the probability density that the next event after H(t ) occurs at time t and has type k , and p(@ ≥ T | H(t )) is the probability that it occurs at some time ≥ T .

In this paper, we will attempt to guess Miss by sampling z values from the posterior distribution

EQUATION

We may obviously drop the second factor of equation ( 3) if (for the given x) it is known to be a constant function of z. In this case, the events are said to be missing at random (MAR). Otherwise, they are missing not at random (MNAR) Little & Rubin (1987) .

It is often intractable to sample exactly from p(z | x), because the first factor p(x z) is complicated (e.g. a neural net). The difficulty is that x and z can be interleaved with each other. As an alternative, we can use normalized importance sampling, drawing many z values from a proposal distribution q(z | x) and weighting them in proportion to p(z|x) q(z|x) . To make it easy to sample from q(z | x), we adopt the following factored definition:

EQUATION

This resembles equation ( 2), but it conditions each proposed unobserved event not only on the history but also on the future F(t i,j ) def = {k i+1 @t i+1 , . . . , k I+1 @t I+1 }. This future consists of all the observed events that happen after t i,j , where in general F(t) def = {k i @t i : t < t i } for any t ∈ [0, T ]. Note the asymmetry with H(t), which includes both observed and unobserved events. q(z | x) can be trained to approximate the target distribution p(z | x), by making

q(• | H, F) ≈ p(• | H, F).

We can sample z from q(z | x) in chronological order: for each 0 ≤ i ≤ I in turn, draw a sequence of J i unobserved events that follow the observed event k i @t i . This sequence ends (thereby determining J i ) if the next proposed event would have fallen after t i+1 and thus is preempted by the observed event k i+1 @t i+1 . Specific p(x z) and q(z | x) distributions will be introduced below.

2.1. THE NEURAL HAWKES PROCESS:

As our generative model of complete event streams Comp in equation ( 2), we need a multivariate point process model. We choose the neural Hawkes process (Mei & Eisner, 2017) , which has proven flexible and effective at modeling many real-world event streams.

Given the history H(t) of all events before time t, the process defines an intensity λ k (t | H(t)) ∈ R ≥0 , which may be thought of as the instantaneous rate of events of type k at time t. More precisely, as dt → 0 + , the number of events of type k occurring in the interval [t, t + dt), divided by dt, approaches λ k (t | H(t)). If no event of any type occurs in this interval (which becomes almost sure as dt → 0 + ), one may still occur in the next interval [t + dt, t + 2dt), and so on.

The intensity functions λ k (t | H(t)) are continuous on intervals during which no event occurs (note that H(t) is constant on such intervals). They jointly determine a distribution over the time of the next event after H(t), as used in every factor of equation (2). As it turns out (Mei & Eisner, 2017) ,

EQUATION

We can therefore train the parameters θ of the λ k functions by maximizing log-likelihood on training data. Each datum is a complete event stream x z over some time interval [0, T ]. In practice, we stop training early when log likelihood stops increasing on held-out development data.

The neural Hawkes process specifically parametrizes

λ k (t | H(t)) as λ k (t | H(t)) = f k (v k h(t)) (6a) f k (x) = s k log(1 + exp(x/s k ))

(6b) The vector h(t) ∈ (-1, 1) D summarizes (t, H(t)). It is the hidden state of a continuous-time LSTM (Mei & Eisner, 2017 ) that read the events in H(t) as they happened, and then waited until time t. The state of such an LSTM evolves endogenously as it waits for the next event, so the timing of the past events has been incorporated into the state h(t).

3. PARTICLE METHODS:

We now describe our particle-based methods for imputing missing z (equation ( 3)). The details are spelled out in Algorithm 1 and Appendix A. For intuition, Figure 1 walks through part of an example.

Algorithm 1 is a Sequential Monte Carlo (SMC) approach. It returns an ensemble of weighted particles Z M = {(z m , w m )} M m=1 . Each particle z m is sampled from the proposal distribution q(z | x), which is defined to support sampling via a sequential procedure that draws one unobserved event at a time. The corresponding w m are importance weights, which are defined as follows, but which are also built up one factor at a time within Algorithm 1:

EQUATION

The numerator of (7b) is given by equations ( 1)-( 2).2 Equation (3) implies that if we could set q(z | x) equal to p(z | x), so that the particles were IID samples from the desired posterior distribution, then all of the wm would be equal and thus w m = 1/m for all m. In practice, q will not equal p, but will be easier than p to sample from. To correct for the mismatch, the importance weights w m are higher for particles that q proposes more rarely than p would.

The distribution formed by the ensemble, p(z), approaches p(z | x) as M → ∞ (Doucet & Johansen, 2009) . Thus, for large M , the ensemble may be used to estimate the expectation of any function f (z), via

EQUATION

In the subsections below, we will describe two specific proposal distributions q that are appropriate for the neural Hawkes process. These distributions define intensity functions λ q over time intervals.

The trickiest part of Algorithm 1 (at line 32) is to sample the next unobserved event from the proposal distribution q. Here we use the thinning algorithm (Lewis & Shedler, 1979; Liniger, 2009; Mei & Eisner, 2017) . Briefly, this is a rejection sampling algorithm whose own proposal distribution uses a constant intensity λ * , making it a homogeneous Poisson process (which is easy to sample from). A event proposed by the Poisson process at time t is accepted with probability λ q (t)/λ * ≤ 1. If it is rejected, we move on to the next event proposed by the Poisson process, continuing until we either accept such an unobserved event or are preempted by the arrival of the next observed event.

3.1. NEURAL HAWKES PARTICLE FILTERING:

We already have a neural Hawkes process p that was trained on complete data. This neurally defines an intensity function λ p k (t | H(t)) for any history H(t) of events before t and each event type k. The simplest proposal distribution uses precisely this process to draw the unobserved events. More precisely, for each i = 0, 1, . . . , I, for each j = 0, 1, 2, . . ., we let the next event k i,j+1 @t i,j+1 be the first event generated by the competing intensity functions λ k (t | H(t)) over the interval t ∈ (t i,j , t i+1 ), where H(t) consists of all observed and unobserved events up through index i, j . If no event is generated on this interval, then the next event is k i+1 @t j+1 . This method is implemented by Algorithm 1 with smooth = false.

3.2. NEURAL HAWKES PARTICLE SMOOTHING:

The neural Hawkes process p(x z) only offers us λ k (t | H(t)). Extra machinery is needed to condition on F(t) as well.

We use a right-to-left continuous-time LSTM, whose details will be shown shortly in Appendix B, to summarize the future F(t) for any time t into another hidden state vector h(t) ∈ R D . Then we parameterize the proposal intensity using an extended variant of equation (6a):

EQUATION

This extra machinery is used by Algorithm 1 when smooth = true. Intuitively, the left-to-right h(t), as explained in Mei & Eisner (2017) , is supposed to learn sufficient statistics for predicting the future as it reads the complete history H(t), while the right-to-left h(t) carries back observed information from the future F(t) in order to correct any mistaken posterior beliefs that h(t) carries.

The right-to-left LSTM has the same architecture as the left-to-right LSTM used in the neural Hawkes process (Mei & Eisner, 2017) , but a separate parameter vector. For any time t ∈ (0, T ), it arrives at h(t) by reading only the observed events {k i @t i : t < t i ≤ t I }, i.e., F(t), in reverse chronological order. Formulas are given in Appendix B.

This is very similar to the forward-backward algorithm (Rabiner, 1989) or Kalman smoothing (Rauch et al., 1965) . However, it is a trained approximation rather than an exact method because of the complicated neural models involved; see Lin & Eisner (2018) for careful discussion of the connection. Regardless of the chosen model, particle smoothing is to particle filtering as Kalman smoothing is to Kalman filtering (Kalman, 1960; Kalman & Bucy, 1961) .

3.2.1. LEARNING TO PROPOSE:

Training the proposal distribution q(z | x) means learning its parameters φ, namely the parameters of the right-to-left LSTM together with matrix V. We are interested in minimizing the Kullback-Leibler (KL) divergence between q(z | x) and p(z

EQUATION

) where log q(z | x) is analogous to equation ( 5).The gradient of exclusive KL divergence is:

EQUATION

where p(x z) is given in equation ( 5) and p miss (z | x z) is assumed known to us for any given pair of x and z.

Minimizing inclusive KL divergence aims at high recall-q(z | x) is adjusted to assign high probabilities to all of the good hypotheses (according to p(z | x)). Conversely, minimizing exclusive KL divergence aims at high precision-q(z | x) is adjusted to assign low probabilities to poor reconstructions, so that they will not be proposed. We seek to minimize the linearly combined divergence Div(p q) = βKL(p q) + (1 -β)KL(q p) with β ∈ [0, 1] (12) and our (Adam) training is early-stopped when the divergence stops decreasing on the held-out development set. When tuning our system (Appendix E.2), β = 1 gave the best results (perhaps unsurprisingly).

But how do we measure these divergences between q(z | x) and p(z | x)? Of course, we actually want the expected divergence when the observed sequence x ∼ p. We pretend that the data are distributed identically to the model p, and thus we sample x from our training examples. For the exclusive divergence, we sample z ∼ q(• | x) from our proposal distribution; notice that optimizing q here is essentially the REINFORCE algorithm Williams (1992) . For the inclusive divergence, we again pretend that the data are distributed identically to p, so instead of sampling a value of z ∼ p(• | x), we can simply use the ground-truth z that happened to occur with this x. We obtain each example with its ground truth z by starting with a fully observed sequence Comp and sampling a partition into Obs = x, Miss = z from the known missingness mechanism p miss .

Appendix F discusses situations where training on incomplete data by EM is possible.

4. A LOSS FUNCTION AND DECODING METHOD:

It is sometimes useful to find a single hypothesis ẑ that minimizes the Bayes risk, i.e., the expected distance from the unknown ground truth z * . This procedure is called minimum Bayes risk (MBR) decoding and can be approximated with our ensemble of weighted particles:

EQUATION

where D(z * , z) is the distance between z * and z. We now propose a specific distance function D.

4.1. OPTIMAL TRANSPORT DISTANCE:

The distance between z and z * is defined as the minimum cost to edit z into z * . To accomplish this edit, we must identify the best alignment-a one-to-one partial matching a-of the events in the two sequences. We require any two aligned events to have the same type k. An alignment edge between a predicted event at time t (in z) and a true event at time t * (in z * ) incurs a cost of |t -t * | to move the former to the correct time. Each unaligned event in z incurs a deletion cost of C delete , and each unaligned event in z * incurs an insertion cost of C insert . Thus,

EQUATION

where A(z * , z) is the set of all possible alignments between z * and z. Finding this distance (and its corresponding alignment a) is similar to finding the edit distance or dynamic time warping between two discrete-time sequences, and a similar dynamic programming algorithm is presented in Algorithm 2 of Appendix C. We set insertion and deletion costs to be the same, i.e., C insert = C delete = C, so that the distance is symmetric.

4.2. APPROXIMATE MBR DECODING:

Since aligned events must have the same type, the MBR problem decomposes into separately choosing a set ẑ(k) of type-k events for each k = 1, 2, . . . , K, based on the particles' sets z

(k)

m of type-k events. Thus, we simplify the presentation by omitting (k) throughout this section.

Finding the optimal set ẑ appears to be NP-hard, by analogy with the Steiner string problem. It involves searching over the infinitely many z ∈ Z. Fortunately, the optimal transport distance D define in section 4.1 warrants:

EQUATION

where P(z) is the power set of any given z-the set of subsequences of z (including empty sequence and z itself). That is to say, there exists one subsequence of z that achieves the minimum Bayes risk (or the minimum weighted distance).

Why this is true? Let's suppose t is the time of one event in z * , and it is aligned to one event at t 1 in z . Recall that the alignment cost between these two events is |t -t 1 |, so we can simply change z * by moving t to t 1 to minimize its total alignment cost. Then what if it is aligned to events at t 1 < t 2 ? While t < t 1 , increasing t always decreases its total alignment cost |t -

t 1 | + |t -t 2 | = t 1 + t 2 -2t.

Similarly, we should decrease t if t > t 2 . When t 1 ≤ t ≤ t 2 , the alignment cost is fixed at t 2 -t 1 . Therefore, the total alignment cost of this event is always minimized if t is moved to t 1 or t 2 . This argument easily generalizes to the cases where t is aligned to more aligned events. The generalization to the cases that each alignment has a weight w is also straightforward, and the full proof can be found in Appendix D.1. Now we have simplified this decoding problem as a combinational optimization problem:

EQUATION

which can be approximately solved although it is still NP-hard to exactly solve.

Our algorithm (details in Algorithm 3 of Appendix D) finds ẑ by iteratively (1) finding its optimal alignment a m with each z m (called Align Phase) using this method of section 4.1, and then (2) going through the following phases-each of them will update ẑ and decreases the weighted distance

M m=1 w m D(z m , ẑ, a m ) which is an upper bound of M m=1 w m D(z m , ẑ):

Move Phase Thanks to Theorem 1, given fixed |ẑ| and {a m } M m=1 , M m=1 w m D(z m , ẑ, a m ) can be minimized by simply moving each t i in the current ẑ to one of the times, if there is any, which it is aligned to. Otherwise, we keep t i unchanged.

Delete Phase Then we may delete any event in ẑ and its associated element in each a m if M m=1 w m D(z m , ẑ, a m ) is decreased afterwards. Note that deletion of each event does not depend on one another, so the weighted distance can actually be minimized by independently checking each event in ẑ. The phase tends to discard the events that are aligned to far-apart times (or nowhere). Insert Phase Then we may further decrease M m=1 w m D(z m , ẑ, a m ) by inserting to ẑ some events that can be aligned, at low cost, to those in z m which would have been left alone otherwise. Thanks again to Theorem 1, such fortunate events can be chosen from M m=1 z m .

Note that even though we used D(z m , ẑ, a m ) multiple times in the above description, only the one time in Move Phase needs to actually call the dynamic programming algorithm (section 4.1) to get a m . Details of the full algorithm (and its theoretical guarantee) can be found in Appendix D.

5. EXPERIMENTS:

We compare our particle smoothing method with the particle filtering baseline on multiple real-world and synthetic datasets. See Appendix E for training details (e.g., hyperparameter selection).

5.1. DATASETS:

The datasets that we use in this paper range from short sequences with mean length 15 to long ones with mean length > 300. For each of the datasets, we have fully observed data that can be used for training. For each dev and test example, we held out some events from the fully observed sequence, so we present the x part as input to the proposal distribution but we also know the z part for evaluation purposes. The dataset and preparation details can be found in Appendix E.

Synthetic Datasets:

We first checked that we could successfully impute unobserved events that are generated from known distributions. That is, when the generating distribution actually is a neural Hawkes process, could our method outperform the particle filtering in practice? Is the performance consistent over multiple datasets drawn from different processes? To investigate this, we synthesized 10 datasets, each of which was drawn from a different neural Hawkes process with randomly sampled parameters.

Elevator System Dataset (Crites & Barto, 1996) . A multi-floor building is often equipped with multiple elevator cars that follow cooperative strategies to transport passengers between floors (Lewis, 1991; Bao et al., 1994; Crites & Barto, 1996) . Observing the activities of some of the cars might help us impute those of the others. This domain is particularly interesting because it is representative of many real-world cooperative (or competitive) scenarios.

New York City Taxi Dataset (Whong, 2014) . Each medallion taxi of New York City forms a sequence of time-stamped pick-up and drop-off events in the five boroughs (i.e., Manhattan, Brooklyn, Queens, The Bronx, and Staten Island). Having observed a sequence of drop-off events, it is interesting to see if the proposed method is able to impute the pick-up events (Figure 1 ).

5.2. EVALUATION AND RESULTS:

First, as an internal check, we measure how probable each ground truth reference z * is under the proposal distribution constructed by each method, i.e., log q(z * | x). The results on the 12 datasets are displayed in Figure 2 .

We now make predictions by MBR decoding. Figure 3 plots the improved performance of neural Hawkes particle smoothing (red) vs. particle filtering (blue). 3 It shows the optimal transport distance broken down by (a) how well it is doing at predicting which events happen-measured by the total number of insertions and deletions (x-axis); and (b) how well it is doing at predicting when those events happen-measured by the total move cost (y-axis). Different choices of C yield different trade-offs between these two metrics. Intuitively, when C ≈ 0, the decoder is free to insert and Figure 3 : Optimal transport distance of neural Hawkes particle smoothing ( ) vs. particle filtering ( ) on test data, broken down by the number of insertions and deletions (x-axis) and the total move/alignment cost (yaxis). Both axes are normalized by the number of true events in the test dataset, N n=1 |z * n|. On each dataset, for each C, the achieves an optimal transport distance that is broken down to the x-axis and y-axis. The pointing to indicates the gradient direction along which the D function would like us to improve for this C value. The starting from this shows the actual improvement obtained by switching to particle smoothing (which is, indeed, an improvement because it has positive dot product with

). Overall, the Pareto frontier (convex hull) of the symbols dominates the Pareto frontier of the symbols, lying everywhere to its left. delete event tokens, so none of them is moved. As the cost C increases, fewer event tokens are inserted or deleted (so lower cost on the x-axis), but more event tokens are aligned (so higher cost on the y-axis).

6. DISCUSSION:

Our technical contribution is threefold. First of all, as far as we know, we are the first to develop general sequential Monte Carlo methods (Moral, 1997; Liu & Chen, 1998; Doucet et al., 2000; Doucet & Johansen, 2009) to approximate a target posterior distribution q(z | x) that is based on a neural model of the complete sequence x z (such as the neural point processes of Du et al. (2016) and Mei & Eisner (2017) ). Most similar to our work is Linderman et al. (2017) 's sequen-tial Monte Carlo method. They modeled complete sequences by a Hawkes process with latent variables-substituting our neural Hawkes process would obtain exactly our particle filtering method (section 3.1). However, our full method has a particle smoother that takes future observations into account while proposing events at any time t, which is not considered in Linderman et al. (2017) . Shelton et al. (2018) developed a reversible jump Markov chain Monte Carlo (MCMC) sampler (Green & Hastie, 2009) for a Hawkes process with latent variables and it allows future observations to "reach back and suggest possible events earlier in the timeline". But their method takes advantage of the Poisson cluster process representation of Hawkes process-each past event generates a Poisson process with exponentially decaying intensity and each new event belongs to one of the processes. Such a representation cannot be established for a neural point process under which the sequence of all past events as a whole generates only one Poisson process, so their method cannot adapt to the complicated neural model that we consider in this work. Other work that infers unobserved events in continuous time also assumes that complete sequences follow a model that is easier to handle than a neural model, including those based on Markov jump processes (Rao & Teh, 2012; 2013) and continuous-time Bayesian networks (Fan et al., 2010) . Lin & Eisner (2018) design a neural particle smoothing algorithm that performs dynamic-programming-style approximate inference on a neural sequential model. But they only work on discrete-time sequences; our neural Hawkes particle smoothing can be seen as a continuous-time generalization of their method.

Secondly, we define the optimal transport distance between event sequences, which is a provably valid metric. It is a Wasserstein distance (Villani, 2008) , or Earth Mover's distance (Kantorovitch, 1958; Levina & Bickel, 2001) , generalized to unnormalized "distributions". There is more than one way to make this generalization, and this is still a subject of active research (Benamou, 2003; Chizat et al., 2015; Frogner et al., 2015; Chizat et al., 2018) . Our definition allows event insertion and deletion while aligning them, but these operations can only apply to an entire event-we cannot align half of an event and delete the other half. Due to these constraints, a dynamic programming rather than linear programming (relaxation) is needed to find the optimal transport. Xiao et al. (2017) also proposed an optimal transport distance between event sequences and it also allows event insertion and deletion. However, their insertion and deletion cost depends on where it is on the time axis while ours doesn't. 45 Our optimal transport distance is similar to the dynamic time warping (DTW) and its extensions (Sakoe & Chiba, 1971; Listgarten et al., 2005) ,which is also a method to calculate an optimal match between sequences. But major difference indeed exists: 1) DTW aligns each event to at least one event and does not allow insertion or deletion, while ours aligns each to at most one and allows insertion and deletion; 2) DTW forces the first-to-first and last-to-last alignment while ours does not; 3) DTW does not allow crossing edges while our cost function is such that within the matching for type k, there is always an optimal solution with no crossing edges (although another equally good solution may have some).

Last but not least, we design an algorithm to find a sequence of unobserved events whose expected distance to an unknown ground truth reference is approximately minimized, and this reference follows a distribution approximated by a set of (weighted) particles. This problem is similar to finding a consensus representation (Steiner sequence) of a set of sequences, which is described in Gusfield (1997) through the concept of multiple sequence alignment (MSA) (Mount, 2004 ) in computational biology. The latter is usually solved by progressive alignment construction using a guide tree (Feng & Doolittle, 1987; Larkin et al., 2007; Notredame et al., 2000) and iterative realignment of the initial sequences with addition of new sequences to the growing MSA (Hirosawa et al., 1995; Gotoh, 1996) . However, these methods cannot be directly applied to our setting because each event in our sequence is a k i @t i pair, not just a discrete type k i -when two events are aligned, we would like their times as well as their types to match.

On multiple synthetic and real-world datasets, our method turns out to be effective at inferring the ground-truth sequence of unobserved events. The improvement of particle smoothing upon particle filtering is substantial and consistent, showing the benefit of training a proposal distribution.

Appendices:



A SEQUENTIAL MONTE CARLO DETAILS:

Our main algorithm is presented as Algorithm 1. It covers both particle filtering and particle smoothing, with optional multinomial resampling.

In this section, we also provide some further notes that are not covered in the pseudocode.

Managing LSTM state information In Algorithm 1, when we push events to stacks H m and F, we update the respective LSTM's configurations (including gates, cell memories and states), and we revert these updates when we pop events from F. These operations facilitate the computation of intensities λ p k (t) and λ q k (t).

Integral computation Computing the weight w m requires handling the integral on line 28 and line 40 of Algorithm 1. We use the same trick in section B.2 of Mei & Eisner (2017) that gives an unbiased estimate of the integrals by evaluating (t i,j -t i,j-1 )λ p (t) and (t i,j -t i,j-1 )λ q (t) at a random t ∼ Unif(t i,j-1 , t i,j ). We can draw N samples instead of one, and this Monte Carlo algorithm will average over these samples to reduce the variance of this noisy estimator.

Choice of λ * How do we construct the upper bound λ * (line 33 of Algorithm 1)? For particle filtering, we follow the recipe in B.3 of Mei & Eisner (2017) : we can express

λ * = f k (max t g 1 (t) + . . . + max t g n (t)) where each summand v kd h d (t) = v kd • o id • (2σ(2c d (t)) -1) is upper-bounded by max c∈{c id ,c id } v kd • o id • (2σ(2c) -1)

. Note that the coefficients w kd may be either positive or negative.

For particle smoothing, we simply have more summands inside f k so λ * = f k (max t g 1 (t) + . . . + max t g n (t) + max t ḡ1 (t) + . . . + max t ḡn (t)) where each extra summand

u kd hd (t) = u kd • o id • (2σ(2c d (t)) -1) is upper-bounded by max c∈{c id ,c id } u kd • o id • (2σ(2c) -1)

and each u kd is the d-th element of vector v k V (equation ( 9)). Note that the o id , c id , cid of newly added summands ḡ are actually from the right-to-left LSTM while those of g are from the left-to-right LSTM. We only use the same notation here for presentation simplicity.

Missing data factors in p Recall that the joint model (1) includes a factor p miss (Miss = z | Comp = x z), which appears in the numerator of the unnormalized importance weight (7b). Regardless of the form of this factor, it could be multiplied into the particle's weight wm at the end of sampling (line 15). However, Algorithm 1 assumes a missingness mechanism where the missingness of each event k@t depends only on that event and preceding events,6 so that p miss (Miss = z | Comp = x z) factors as

∈indices(x z) p miss ((k @t ∈ Miss) = (k @t ∈ z) | {k @t : ≤ }) (17)

Algorithm 1 can thus incrementally incorporate the subfactors of equation ( 17), and does so at line 30.

Optional missing data factors in q We can optionally improve the particle filtering proposal intensities to incorporate the p miss factor discussed above (in which case that factor will be multiplied into the denominator of (7b) and not just the numerator). This makes q(z | x) better match p(z | x): it means we will rarely posit an unobserved event that would rarely have gone missing.

Specifically, if a completed-data event k@t would have probability r k (t | H(t)) of going missing given the preceding events H(t), it is wise to define

λ q k (t | H(t)) = λ p k (t | H(t)) • r k (t | H(t)

). We include this extra r k factor in our experiments (section 5), although it is not shown in Algorithm 1. It is particularly simple in our experiments, where r k is constant at 1 or 0 depending on k. In other words, some event types k are deterministically missing, while others are never missing and thus we never propose them as part of z.

The above considers the particle filtering case. In the case of particle smoothing, we have already tried to ensure by other means that the proposal distribution will incorporate p miss . That is because section 3.2.1 aims to train λ q k (t | H(t), F(t)) so that the resulting q(z | x) ≈ p(z | x), and the posterior distribution p(z | x) does condition on the missingness of z. However, if the r k factor is known, why not include it explicitly instead of having to train the BiLSTM to mimic it? Thus it can be convenient to modify the right-hand side of equation ( 9) to include a factor of r k . This yields a more expressive and better-factored family of proposal distributions: missingness is now handled by the known r k factor and the BiLSTM does not have to explain it.

Modifying equation ( 9) in this way is particularly useful in the special case r k = 0 (i.e., event type k is never missing and should not be proposed). There, it enforces the hard constraint that λ q k = 0 (something that the BiLSTM by itself could not achieve); and since this constraint is enforced regardless of the BiLSTM parameters, the events of type k appropriately become irrelevant to the training of the BiLSTM, which can focus on predicting other event types. We do this in our experiments.

B RIGHT-TO-LEFT CONTINUOUS-TIME LSTM:

Here we give details of the right-to-left LSTM from section 3.2. At each time t ∈ (0, T ), its hidden state h(t) is continually obtained from the memory cells c(t) as the cells decay:

EQUATION

) where the interval (t i-1 , t i ) has consecutive observations k i-1 @t i-1 and k i @t i as endpoints. At t i , the continuous-time LSTM reads k i @t i and updates the current (decayed) hidden cells c(t) to new initial values c i-1 , based on the current (decayed) hidden state h(t i ), as follows:

7 i i-1 ← σ (W i k i + U i h(t i ) + d i ) (19a) f i-1 ← σ (W f k i + U f h(t i ) + d f ) (19b) z i-1 ← 2σ (W z k i + U z h(t i ) + d z ) -1 (19c) o i-1 ← σ (W o k i + U o h(t i ) + d o ) (19d) c i-1 ← f i-1 c(t i ) + i i-1 z i-1 (20a) ci-1 ← f i-1 ci + īi-1 z i-1 (20b) δ i-1 ← f (W d k i + U d h(t i ) + d d ) (20c)

The vector k i ∈ {0, 1} K is the i th input: a one-hot encoding of the new event k i , with non-zero value only at the entry indexed by k i . Then, c(t) is given by ( 21), which continues to control h(t) except that i has now decreased by 1). c(t)

EQUATION

On the interval [t i-1 , t i ), c(t) follows an exponential curve that begins at c i-1 (in the sense that lim t→t - i c(t) = c i-1 ) and decays, as time t decreases, toward ci-1 .

C OPTIMAL TRANSPORT DISTANCE DETAILS:

In this section, we first present the detailed algorithm of finding the OTD and its corresponding alignment, and then prove it is a valid metric.

C.1 ALGORITHM DETAILS:

The details of how to find optimal transport distance is presented in Algorithm 2.

C.2 PROOF THAT OPTIMAL TRANSPORTATION DISTANCE IS A VALID METRIC:

We have defined that the optimal transportation distance (OTD), and we prove that it is a valid metric here.

Algorithm 1 Sequential Monte Carlo -Neural Hawkes Particle Filtering/Smoothing:

Input: observed seq. x = k 0 @t 0 , . . . , k I+1 @t I+1 with k 0 = BOS, t 0 = 0, k I+1 = EOS, t I+1 = T ; model p; missingness mechanism p miss ; proposal distribution q; number of particles M ; boolean flags smooth and resample Output: collection {(z 1 , w 1 ), . . . , (z M , w M )} of weighted particles 1: procedure SEQUENTIALMONTECARLO(x, p, p miss , q, M, smooth, resample)

2:

for m = 1 to M : for i = I downto 0 : stack of all future observed events 7: push k i+1 @t i+1 onto F as we reach these events, we'll pop from F and push onto Hm (∀m) 8: else use particle filtering instead 9:

F ← ignored special value if we're not using the future stack; unaffected by pop operation 10:

for i = 0 to I : observe present event ki@ti, then propose unobserved events on interval (ti, ti+1) 11:

for m = 1 to M : for m = 1 to M : often draws multiple copies of good (high-weight) particles, 0 copies of bad ones 18:

zm ∼ Categorical({z m → w m / M m=1 w m } M m=1 )

19:

for m = 1 to M : 

23::

q gives info to define function λ q k (t)

def = λ k (t | Hm, F), or λ q k (t) = λ p k (t) if F = ignored 24:

these functions consult state of left-to-right LSTM that's read Hm & right-to-left LSTM that's read F 25:

we also define the total intensity functions λ p (t)

def = K k=1 λ p k (t) and λ q (t) def = K k=1 λ q k (t) 26: j ← 0; k i,j ← k i ; t i,j ← t i

ready to observe i th event of x 27:

while true : one iteration adds event at index i, j (for j = 0, 1, 2, . . . until we break out of loop) 28:

w m ← w m • λ p ki,j (t i,j ) exp (- ti,j t =ti,j-1 λ p (t )dt )

new factor in numerator p of (7b) 29: push k i,j @t i,j onto H m ; t ← t i,j event just generated by p now becomes part of history 30:

w m ← w m • p miss ((k i,j @t i,j ∈ Miss) = (j > 0) | H m )

new factor in numerator p of (7b) 31:

Now draw possible missing event at index i, j +1 ; we'll loop back and add it if it falls in (ti, ti+1) 32: repeat thinning algorithm (see Mei & Eisner, 2017) 33: find any λ * ≥ sup {λ q (t ) : t ∈ (t, t i+1 )} e.g., old λ * still works if i unchanged 34:

draw ∆ ∼ Exp(λ * ), u ∼ Unif(0, 1) 35: t += ∆ time of next proposed event (before thinning) 36: if t ≥ t i+1 :

proposed event falls outside (ti, ti+1), where ki+1@ti+1 is top element of F 37: return done adding events on [ti, ti+1); time to break out of while loop 38:

until uλ * ≤ λ q (t) thinning: accept proposal with prob λ q (t) λ * ≤ 1 39:

j ← j + 1; t i,j ← t; k i,j ∼ Categorical({k → λ q k (t) λ q (t) } K k=1 ); append k i,j @t i,j to z m 40:

w m ← w m / λ q ki,j (t i,j ) exp(-ti,j t =ti,j-1 λ q (t )dt )

new factor in denominator q of (7b)

It's trivial that OTD is non-negative, since movement, deletion and insertion costs are all positive.

It's also trivial to prove that the following statement is true:

EQUATION

where z 1 and z 2 are two sequences. If z 1 is not identical to z 2 , the distance of them must be larger than 0 since we have to do some movement, insertion or deletion to make them exactly matched, so the right direction of equation ( 22) holds. If the distance between z 1 and z 2 is zero, which means for k ← 1 to K : 4: for i * ← 1 to I * : to match reference of length i * 17: if P î,i * = insert :

d (k) , a (k) ← DYNAMICPROGRAMMING(ẑ (k) , z * (k) ) 5: d ← d + d (k) ; a ← a ∪ a (k) 6: return d, a 7: procedure DYNAMICPROGRAMMING(z * (k) , ẑ(k) ) 8: Î ← |ẑ (k) |; I * ← |z * (k) | ẑ(k) =

D insert ← D î-1,i * + C insert if an event token at tî is inserted to z * (k) 18: D delete ← D î,i * -1 + C delete if the event token at t * i * is deleted from z * (k) 19: D move ← D î-1,i * -1 + | tî -t * i * | if the event at t * i * of z * (k) is aligned to event at tî of ẑ(k) 20: D î,i * ← min{D insert ,

token t * i * is deleted. 25: î ← î -1 26: if P î,i * = delete : a token at tî is inserted 27: i * ← i * -1 28: if P î,i * = move : token t * i * is aligned to tî 29: î ← î -1; i * ← i * -1 30: a ← a ∪ {( tî , t * i * )} 31:

return D Î,I * , a they are already matched without any operations, z 1 and z 2 must be identical, thus the left direction of equation ( 22) holds.

OTD is symmetric, that is,

D(z 1 , z 2 ) = D(z 2 , z 1 ), if we set C insert = C delete .

Suppose that a is an alignment between z 1 and z 2 . It's easy to see that the only difference between D(z 1 , z 2 , a) and D(z 2 , z 1 , a)8 is that the insertion and deletion operations are exchanged. For example, if we delete a token t i ∈ z 1 when calculating D(z 1 , z 2 , a), we should insert a token at t i to z 2 when calculating

EQUATION

)

Finally let's prove that OTD satisfies triangle inequality, that is:

EQUATION

where z 1 , z 2 and z 3 are three sequences. This property could be proved intuitively. Suppose that the operations on z 1 with minimal costs to make z 1 matched to z 2 are denoted by o 1 , o 2 , . . . , o n1 , and those on z 2 to make z 2 matched to z 3 are denoted by o 1 , o 2 , . . . , o n2 . o i could be a deletion, insertion or movement on a token. To make z 1 matched to z 3 , one possible way, which is not necessarily the optimal, is to do o 1 , o 2 , . . . , o n1 , o 1 , o 2 , . . . , o n2 on z 1 . Since the total cost is the accumulation of the cost of each operation, and the operations on z 1 above to make z 1 matched to z 3 might not be optimal, the triangle inequality equation ( 25) holds.

D APPROXIMATE MBR DETAILS:

In this section, we first prove the Theorem 1 in section 4.2, and then show the detailed algorithm to find the decode.

D.1 THEOREM PROOF AND RELATED:

We have a claim in section 4.2 that:

EQUATION

where P(z) is the power set of any given z-the set of subsequences of z (including empty sequence and z itself). That is to say, there exists one subsequence of z that achieves the minimum Bayes risk (or the minimum weighted optimal transport distance).

and we prove it in this section:

Proof. Here we assume that there is only one type of event. Since the distances of different types of events are calculated separately, our conclusion is easy to be extended to the general case.

Suppose z * is an optimal decode, that is, If z * ∈ P(z ), the proof is done. If not, suppose there exists a token at t i / ∈ z . We use t l ∈ z (t r ∈ z ) to denote the token in z that is left (right) and nearest to t i . 9 We will show that if we move t i around, as long as t i ∈ [t l , t r ], the weighted optimal transport distance, i.e. Let's use r(t) to indicate the weighted transport distance of z * with fixed alignment if we move t i to t, that is,

r(t) def = M m=1 w m D(z m , z * (t), a * m ),

where z * (t) is the sequence z * with t i moved to t. Because z * (t i ) is an optimal decode, and a * m is the optimal alignment for z * (t i ), we should have

r(t i ) = min t r(t).

Note that the transport distance is comprised of three parts: deletion, insertion and alignment costs. Since every a * m is fixed, if we change t, only the alignment cost that related to token t will affect r(t). This part of r(t) is linear to t, since we have a constraint t ∈ [t l , t r ], which guarantees that it will not cross any other tokens in z .

Since r(t) is linear to t ∈ [r l , t r ] and r(t) gets minimized at t i ∈ (t l , t r ), we could conclude that Since r(t) is the upper bound of the weighted optimal transport distance M m=1 w m D(z m , z * (t)), which also gets the same minimal value at t i ∈ (t l , t r ) as r(t), we could conclude that

r(t) = r(t i ) = Const, ∀t ∈ [t l , t r ].

M m=1 w m D(z m , z * (t)) = M m=1 w m D(z m , z * (t i )) = Const, ∀t ∈ [t l , t r ].

Therefore we could move token t i to either t l or t r without increasing the Bayes risk. We could do this movement for each t i / ∈ z to get a new decode ẑ ∈ P(z ), which is also an optimal decode.

D.2 ALGORITHM DETAILS:

The detailed algorithm is presented in Algorithm 3.

E EXPERIMENTAL DETAILS:

In this section, we elaborate on the details of data generation, processing, and experimental results.

In all of our experiments, the distribution p is trained on the complete (uncensored) version of the training data. The system is then asked to complete the incomplete (censored) version of the test (or dev) data. For particle smoothing, the proposal distribution is trained using both the complete and incomplete versions of the training data, as explained at the end of section 3.2.1.

In each experiment, the missingness mechanism is defined to be as

EQUATION

meaning that each event in the complete stream x z is independently censored with probability ρ k that only depends on its event type k. Our experiments on imputing missing data assume that both the missingness mechanism equation ( 27) and its parameter vector ρ are known, although Appendix F discusses how ρ could be imputed when complete and incomplete data are both available..

For our main experiments in section 5.2, ρ k is either 0 or 1, which means some event types are always observed, while others are always missing. Though this missingness mechanism is deterministic, it can still be regarded as MNAR, because the factor p miss (Miss = z | Comp = x z) in equation ( 3) is not constant. The factor is 1 if z consists of precisely the events in x z that ought to go missing, and 0 otherwise. In other words, p miss here ensures only that z should consist only of events with ρ k > 0. Recall from Appendix A that p miss is also considered by the proposal distribution.

We also experiment with a more typical MNAR setting in Appendix G below, in which events are stochastically missing rather than deterministically missing.

E.1 DATASET STATISTICS:

Table 1 shows statistics about each dataset that we use in this paper.

E.2 TRAINING DETAILS:

We used single-layer LSTMs (Hochreiter & Schmidhuber, 1997) , selected the number D of hidden nodes of the left-to-right LSTM, and then D of the right-to-left one from a small set For the chosen (D, D ) pair on each dataset, we selected β based on the performance on the dev set, and β = 1.0 yields the best performance across all the datasets we use. For learning, we used Adam with its default settings (Kingma & Ba, 2015) .

E.3 SYNTHETIC DATASETS DETAILS:

Each of the ten neural Hawkes processes has its parameters sampled from Unif[-1.0, 1.0]. Then a set of event sequences is drawn from each of them via the plain vanilla thinning algorithm (Mei & Eisner, 2017) . For each of the ten synthetic datasets, we took K = 4 as the number of event types. To draw each event sequence, we first chose the sequence length I (number of event tokens) uniformly from {11, 12, . . . , 20} and then used the thinning algorithm to sample the first I events over the interval [0, ∞). For subsequent training or testing, we treated this sequence (appropriately) as the complete set of events observed on the interval [0, T ] where T = t I , the time of the last generated event.

We generate 5000, 500 and 500 sequences for each training, dev, and test set respectively. For the missingness mechanism, we censor all events of type 3 and 4. In other words, we set ρ k = 0 for k = 1, 2 and ρ k = 1 for k = 3, 4.

E.4 ELEVATOR SYSTEM DATASET DETAILS:

We examined our method in a simulated 5-floor building with 2 elevator cars. During a typical afternoon down-peak rush hour (when passengers go from floor-2,3,4,5 down to the lobby), elevator cars travel to each floor and pick up passengers that have (stochastically) arrived there according to a traffic profile (Bao et al., 1994) . Each car will also avoid floors that already are or will soon be taken care of by the other. Having observed when and where car-1 has stopped (to pick up or drop off passengers) over this hour, we are interested in when and where car-2 has stopped during the same time period. In this dataset, each event type is a tuple of (car number, floor number) so there are K = 10 in total in this simulated 5-floor building with 2 elevator cars.

Passenger arrivals at each floor are assumed to follow a inhomogeneous Poisson process, with arrival rates that vary during the course of the day. The simulations we use follows a human-recorded traffic profile (Bao et al., 1994 ) which dictates arrival rates for every 5-minute interval during a typical afternoon down-peak rush hour. Table 2 shows the mean number of passengers (who are going to the lobby) arriving at floor-2,3,4,5 during each 5-minute interval.

We simulated the elevator behavior following a naive baseline strategy documented in Crites & Barto (1996) . 10 In details, each car has a small set of primitive actions. If it is stopped at a floor, it must either "move up" or "move down". If it is in motion between floors, it must either "stop at the next floor" or "continue past the next floor". Due to passenger expectations, there are two constraints on these actions: a car cannot pass a floor if a passenger wants to get off there and cannot turn until it has serviced all the car buttons in its current direction. Three additional action constraints were made in an attempt to build in some primitive prior knowledge: 1) a car cannot stop at a floor unless someone wants to get on or off there; 2) it cannot stop to pick up passengers at a floor if another car is already stopped there; 3) given a choice between moving up and down, it should prefer moving up (since the down-peak traffic tends to push the cars toward the bottom of the building). Because of this last constraint, the only real choices left to each car are the stop and continue actions, and the baseline strategy always chooses to continue. The actions of the elevator cars are executed asynchronously since they may take different amounts of time to complete.

We repeated the (one-hour) simulation 700 times to collect the event sequences, each of which has around 300 time-stamped records of which car stops at which floor. We randomly sampled disjoint train, dev and test sets with 500, 100 and 100 sequences respectively.

We set ρ k = 0 for k = 1, . . . , 5 and ρ k = 1 for k = 6, . . . , 10, meaning that the events (of arriving at floor 1, 2, . . . , 5) of car 1 are all observed, but those of car 2 are not. The New York City Taxi dataset (section 5.1) includes 189,550 taxi pick-up and drop-off records in the city of New York in 2013. Each record has its medallion ID, driver license and time stamp. Each combination of medallion ID and driver license naturally forms a sequence of time-stamped pick-up and drop-off events. Following the processing recipe of previous work (Du et al., 2016) , we construct shorter sequences by breaking each long sequence wherever the temporal gap between a drop-off event and its following pick-up event is larger than six hours. Since the schedule of different drivers might be very different, it's hard to set a natural t 0 , that is, the time stamp of BOS.

We randomly sampled a month from 2013 and then randomly sampled disjoint train, dev and test sets with 5000, 500 and 500 sequences respectively from that month.

In this dataset, each event type is a tuple of (borough, action) where the action can be either pick-up or drop-off, so there are K = 5 × 2 = 10 event types in total.

We set ρ 1 = 0 and ρ 2 = 1, which means that all drop-off events but no pick-up events are observed.

F MONTE CARLO EM:

We normally assume (section 3.2.1) that some complete sequences are available for training the neural Hawkes process models. If incomplete sequences are also available, our particle smoothing method can be used to (approximately) impute the missing events, which yields additional complete sequences for training. Indeed, if we are willing to make an MAR assumption (Little & Rubin, 1987) ,11 then we can do imputation without modeling the missingness mechanism. Training on such imputed sequences is an instance of Monte Carlo expectation-maximization (MCEM) (Dempster et al., 1977; Wei & Tanner, 1990; McLachlan & Krishnan, 2007) , with particle smoothing as the Monte Carlo E-step, and makes it possible to train with incomplete data only. In the more general MNAR scenario, we can extend the E-step to consider the not-at-random missingness mechanism (see equation (7b) below), but then we need both complete and incomplete sequences at training time in order to fit the parameters of the missingness mechanism (unless these parameters are already known) jointly with those of the neural Hawkes process. Although training with incomplete data is out of the scope of our experiments, we describe the methods and provide MCEM pseudocode in Appendix F.

In this case, we would like to know the probability of the observed data under the target distribution:

EQUATION

where the marginal target distribution p(Obs = x) is abbreviated as p(x). If we propose z from q(z | x), then it can be rewritten as:

EQUATION

Given a finite number M of proposed particles {z m } M m=1 , this expectation can be estimated with empirical average: 

G MORE EXPERIMENTS ON MNAR DATA:

Our experimental setting in section 5.1 and Appendix E is a special case of MNAR where some event types are deterministically censored. In this section, we consider a more typical MNAR setting, where (in the notation of Appendix E) we set ρ k = ρ = 0.5 regardless of the event type k. Then equation ( 27) can be written as p miss (Miss = z | Comp = x z) = (0.5)

|x|+|z| , whose value depends on the length of |z|. This is obviously non-constant over z.

The proposal distribution in this case will be more conservative about proposing missing events,because having a lot of missing events is a posteriori improbable. In other words, p miss as given above falls off with the number of missing events |z|.

We conducted experiments on the real-world datasets in this setting. Again, the method works, with very similar qualitative behavior to before. As shown in Figures 4 and 5 , the particle smoothing method outperforms the filtering baseline (i.e. the neural version of Linderman et al. (2017) ). Note that these figures look very similar to Figures 2 and 3 . This is not too surprising because these two settings both have about half their events censored-thus being about equally difficult to reconstruct. 

:

Our conventions of mathematical notation mainly follow those given byMei & Eisner (2017, section

2).

As discussed earlier, if we are willing to make a MAR assumption, we may omit the pmiss factor of equation (1) because it is constant (though unknown).

We show the 2 real datasets only. The figures for the 10 synthetic datasets are boringly similar to these.

This dependence actually makes the distance unintuitive. For example, on interval [0, 100), the distance between z1 = 1, 3, 4,

and z2 = 3, 4, 5 is |1 -3| + |3 -4| + |4 -5| + |5 -100| = 99 under their definition while it can be small under our definition (e.g. 0.5 if our insertion and deletion cost C = 0.5). The latter seems more natural because these two sequences are almost identical.5 Besides optimal transport, Stein's method(Stein et al., 1972) is also used to find (or bound) the distance of point processes(Schuhmacher & Xia, 2008;Decreusefond et al., 2016).

This assumption could trivially be relaxed to allow it to also depend on the missingness of the preceding events, and/or on the future observed events F(t).

The upright-font subscripts i, f, z and o are not variables, but constant labels that distinguish different W, U and d tensors. The f and ī in equation (20b) are defined analogously to f and i but with different weights.

We abuse the notation a, which we think could represent both the movement from z1 to z2 and from z2 to z1.

We assume that ti is in between two tokens in z , and our proof could be easily extended to the case for which this condition is not satisfied.

We rebuilt the system in Python following the original Fortran code ofCrites & Barto (1996).

It is "almost impossible" to determine from the data whether the MAR assumption holds(Mohan & Pearl, 2018).

Figures and Tables:

Figure 2: Scatterplots of neural Hawkes particle smoothing (y-axis) vs. particle filtering (x-axis). Each point represents a single test sequence, and compares the values of log q(z * | x)/|z * | (i.e., nats per unobserved event)under the two proposal distributions. Larger numbers mean that the proposal distribution is better at proposing the ground truth. Each dataset's scatterplot is converted to a cloud using kernel density estimation, with the centroid denoted by a black dot. A double-arrowed line indicates the improvement of particle smoothing over filtering. For the synthetic datasets, we draw ten clouds on the same figure and show the line for the set where smoothing improves the most. Particle smoothing performs well even on the datasets where particle filtering performs badly. As we can see, the density function is always well concentrated above y = x. That is, this is not merely an average improvement: nearly every ground truth sequence increases in proposal probability!

init weighted particles (zm, wm). History Hm combines zm with a prefix of x 3: z m ← empty sequence; w m ← 1; H m ← empty stack 4:

m) destructively extend zm, wm, Hm with events on [ti, ti+1) 13: pop F pop ki+1@ti+1 from F: it's no longer in the future but in the present 14: if resample : RESAMPLE() optional multinomial resampling replaces all weighted particles 15: return {(z m , w m / M m=1 w m )} M m=1 M particles with weights normalized as in equation (7a) 16: procedure RESAMPLE has access to global variables 17:

zm ; w m ← 1 update particles and their weights 21: procedure DRAWSEGMENT(i, m) has access to global variables 22: p gives info to define function λ p k (t) def = λ k (t | Hm)

Dynamic Programming Algorithm to Find Optimal Transport Distance Input: proposal ẑ; reference z * Output: optimal transport distance d; alignment a 1: procedure OTD(z * , ẑ) 2: d ← 0; a ← empty collection {} 3:

t1, . . . , tÎ and z * (k) = t * 1 , . . . , t * I * 9: D ← zero matrix with (M + 1) rows and (N + 1) columns 10: P ← empty matrix with M rows and N columns back pointers 11: for î ← 1 to Î : transport reference of length 0 to proposal of length î 12: D î,0 ← D î-1,0 + C insert insert t * î = tî to reference (and their prefixes are matched) 13: for i * ← 1 to I * : transport preference of length i * to proposal of length 0 14: D 0,j * ← D 0,j * + C delete delete t * i * (and prefixes are matched) 15:for î ← 1 to Î :proposal prefix of length î 16:

D delete , D move } choose the edit that yields the shortest distance 21: P î,i * ← arg min e∈{insert,delete,move} D e e represents a kind of edition 22: î ← Î; i * ← I * ; a ← empty collection{} 23: while î > 0 and i * > 0 : back trace 24:

m=1 w m D(z m , z * ), will neither increase nor decrease. Suppose a * m = arg min am∈A(zm,z * ) M m=1 w m D(z m , z * , a m ).

Figure 4: Scatterplots with ρ = 0.5. Same comparison as Figure 2.

Find an approximately minimum-Bayes-risk sequence of events Input: collection of weighted particles Z M = {(z m , w m )} DECODEK(Z M )

Competitive games. In competitive games such as StarCraft and poker, a player does not have full information about what her opponents have done (e.g., build mines and train soldiers) or acquired (e.g., certain cards). Correctly imputing what has happened "what I did"

Table content: <table/>

Statistics of each dataset. We write "≈ N " to indicate that N is the average value over multiple datasets of one kind (synthetic); the variance is small in each such case.

Table content: <table><tr><td>DATASET</td><td>K</td><td colspan="2"># OF EVENT TOKENS</td><td/><td colspan="3">SEQUENCE LENGTH</td></tr><tr><td/><td/><td>TRAIN</td><td>DEV</td><td>TEST</td><td>MIN</td><td>MEAN</td><td>MAX</td></tr><tr><td>SYNTHETIC</td><td>4</td><td>≈ 74967</td><td>≈ 7513</td><td>≈ 7507</td><td>10</td><td>≈ 15</td><td>20</td></tr><tr><td>NYCTAXI</td><td>10</td><td>157916</td><td>15826</td><td>15808</td><td>22</td><td>32</td><td>38</td></tr><tr><td>ELEVATOR</td><td>10</td><td>313043</td><td>31304</td><td>31206</td><td>235</td><td>313</td><td>370</td></tr></table>

{16, 32, 64, 128, 256, 512, 1024}  based on the performance on the dev set of each dataset. The best-performing (D, D ) pairs are (256, 128) on Synthetic, (256, 256) on Elevator (256, 256) on NYC Taxi, but we empirically found that the model performance is robust to these hyperparameters.

Table content: <table/>

The Down-Peak Traffic Profile E.5 NEW YORK CITY TAXI DATASET DETAILS

Table content: <table><tr><td>START TIME (MIN)</td><td>00</td><td>05</td><td>10</td><td>15</td><td>20</td><td>25</td><td>30</td><td>35</td><td>40</td><td>45</td><td>50</td><td>55</td></tr><tr><td>MEAN # PASSENGER</td><td>1</td><td>2</td><td>4</td><td>4</td><td>18</td><td>12</td><td>8</td><td>7</td><td>18</td><td>5</td><td>3</td><td>2</td></tr></table>

Peer review (iclr2019):
 The authors propose a particle smoothing approach with an approximate minimum Bayes risk decoder to impute missing events in the Neural Hawkes Process (NHP). The main goal is to address the missing events problem in continuous-time event analysis, which is an important problem in practice. The core idea is within the framework of particle smoothing. To formulate the posterior distribution of the missing event, the authors consider both the left-to-right past events and the right-to-left future events. The paper first applies the NHP to capture both the observed and inferred missing events to learn a representation of the past events, and then uses a similar NHP to learn the representation of the observed events from the future. Based on the two representations, it then formulates the intensity function of the missing events and uses the thinning algorithm to sample different particles. Based on the proposed distribution, the paper also considers to decode a single prediction achieving the Minimum Bayes Risk. Experiments on synthetic datasets with 10 different initializations and two real datasets show that the proposed smoothing approach is better than the filtering baseline. In general, this paper considers an important problem which is under active research in literature recently. However, there are a few weaknesses of the paper that should be addressed. 1. The proposed technique is tightly connected to NHP, which could limit the applicability of the approach to other temporal point processes. The essential idea is similar to Bi-LSTM to learn the representation from both ends of a sequence of asynchronous temporal events. There are several different ways to represent the inter-event time to feed into the network other than NHP. Can the proposed method also be applied to other processes? 2. Within the particle filtering framework, each particle (hypothesis) is weighted by the likelihood of the sequence of observed events under that hypothesis. It turns out that the integral part of Equation 5 does not have an obvious analytical solution under NHP. Then, we first need a set of samples to approximate the likelihood evaluation. Later, we also need to sample particles. I am not quite convinced the computational efficiency of this approach in real applications of practice. Also, there is no analysis either empirically or analytically about the impact of the accumulative sampling errors on the inference performance. Furthermore, to learn the proposed distribution, the paper applies the REINFORCE algorithm under the proposed distribution q. But REINFORCE is known for large variance issue. Given that we already need lots of samples for the likelihood, it is unclear to me how stable the algorithm could be in practice. 3. The experimental evaluation is weak. For particle filtering and smoothing, it is known that the filtering techniques are candidates for solving the smoothing problem but perform poorly when T is large. That*s why it is necessary to develop more sophisticated strategies for good smoothing algorithms. As a result, it is unfair to only compare the smoothing approach with the filtering baseline. Actually, what people really care about is how different techniques can behave in real data to impute realistic missing events. From this perspective, I suggest to use the QQ-plot to evaluate the goodness of fitting on the synthetic dataset. For example, given a sequence of events generated from an independent temporal point process, we can randomly delete events, and then apply different techniques, including Linderman et al. (2017), Shelton et al.(2018), to impute missing events. Finally, we can compare the imputed sequence of events with the groundtruth. In addition, sequential monte carlo approach often suffers from skewed particle issue where one particle gradually dominates all the other particles with no diversity. It is unclear how the proposed approach is able to handle this. One missing related paper is "Learning Hawkes Processes from Short Doubly-Censored Event Sequences" Section 5.2 can be significantly strengthened if comparing with at least one of these approaches. 4. The paper is fairly written. I had some trouble reading back and forth for understanding Figure 1 since it has long caption that is not self-contained. The annotation of Section 2 is also too heavy to quickly skim through to memorize.