Year: 2018

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.

heading: 1 INTRODUCTION
text: 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.
• 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” 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
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 t 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 draws— was 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.

heading: 2 PRELIMINARIES1
text: 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 t Miss , where t 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 = {k0@t0, k1@t1, k2@t2, . . . , kI@tI , kI+1@tI+1}, where each ki ∈ {1, 2, . . . ,K} is an event type and 0 = t0 < t1 < t2 < . . . < tI < tI+1 = T are the times of occurrence in increasing order. This sequence always includes the special boundary events k0@t0 = BOS@0 (“beginning of sequence”) and kI+1@tI+1 = EOS@T (“end of sequence”).
Following each observed event ki@ti for 0 ≤ i ≤ I , we may hypothesize some Ji ≥ 0 unobserved events {ki,1@ti,1, ki,2@ti,2, . . . , ki,Ji@ti,Ji}, where ti < ti,1 < ti,2 < . . . < ti,Ji < ti+1. We may regard the subscript i as a shorthand for (i, 0), so the observed events are x = {ki,j@ti,j : j = 0} and the hypothesized unobserved (missing) events are z = {ki,j@ti,j : j 6= 0}. The hypothesized complete event stream x t 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 < Ji 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 defineH(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 pmiss: p(Obs = x,Miss = z) = p(Comp = x t z) · pmiss(Miss = z | Comp = x t z) (1) Furthermore, our complete data model (such as a neural Hawkes process) will take the factored form
p(Comp = x t z) = ( I∏ i=0 Ji∏ j=0 p(ki,j@ti,j | H(ti,j)) ) · p(@≥T | H(tI,JI )) (2)
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
p(Miss = z | Obs = x) ∝ p(Comp = x t z) · pmiss(Miss = z | Comp = x t z) (3) 1Our conventions of mathematical notation mainly follow those given by Mei & Eisner (2017, section 2).
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(xtz) 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:
q(z | x) = I∏ i=0 (( Ji∏ j=1 q(ki,j@ti,j | H(ti,j),F(ti,j)) ) · q(@≥ ti+1 | H(ti,Ji),F(ti,Ji)) ) (4)
This resembles equation (2), but it conditions each proposed unobserved event not only on the history but also on the futureF(ti,j) def = {ki+1@ti+1, . . . , kI+1@tI+1}. This future consists of all the observed events that happen after ti,j , where in general F(t) def = {ki@ti : t < ti} 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 Ji unobserved events that follow the observed event ki@ti. This sequence ends (thereby determining Ji) if the next proposed event would have fallen after ti+1 and thus is preempted by the observed event ki+1@ti+1.
Specific p(x t z) and q(z | x) distributions will be introduced below.

heading: 2.1 THE NEURAL HAWKES PROCESS
text: 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 afterH(t), as used in every factor of equation (2). As it turns out (Mei & Eisner, 2017),
log p(Comp = x t z) = ∑ `:t`<T log λk`(t` | H(t`))− ∫ T t=0 K∑ k=1 λk(t | H(t))dt (5)
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 t 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)) = fk(v>k h(t)) (6a)
fk(x) = sk log(1 + exp(x/sk)) (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 ZM = {(zm, wm)}Mm=1. Each particle zm 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 wm are importance weights, which are defined as follows, but which are also built up one factor at a time within Algorithm 1:
wm = w̃m∑M m=1 w̃m
(7a) w̃m = p(Obs = x,Miss = zm)
q(zm | x) (7b)
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 w̃m would be equal and thus wm = 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 wm 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
Ep(z|x)[f(z)] ≈ Ep̂[f(z)] = ∑ z p̂(z)f(zm) = M∑ m=1 wmf(zm) (8)
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.

heading: 3.1 NEURAL HAWKES PARTICLE FILTERING
text: We already have a neural Hawkes process p that was trained on complete data. This neurally defines an intensity function λpk(t | H(t)) for any historyH(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 ki,j+1@ti,j+1 be the first event generated by the competing intensity functions λk(t | H(t)) over the interval t ∈ (ti,j , ti+1), whereH(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 ki+1@tj+1. This method is implemented by Algorithm 1 with smooth = false.

heading: 3.2 NEURAL HAWKES PARTICLE SMOOTHING
text: The neural Hawkes process p(x t 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) ∈ RD′ . Then we parameterize the proposal intensity using an extended variant of equation (6a):
λqk(t | H(t),F(t)) = fk(v > k (h(t) + Vh̄(t))) (9)
2As 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).
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 {ki@ti : t < ti ≤ tI}, 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).

heading: 3.2.1 LEARNING TO PROPOSE
text: 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 KullbackLeibler (KL) divergence between q(z | x) and p(z | x). Although p(z | x) is unknown, the gradient of inclusive KL divergence between q(z | x) and p(z | x) is
∇φKL(p‖q) = Ez∼p(z|x)[−∇φ log q(z | x)] (10) where log q(z | x) is analogous to equation (5).The gradient of exclusive KL divergence is:
∇φKL(q‖p) = Ez∼q(z|x)[∇φ ( 1
2 (log q(z | x)− log p(x t z)− pmiss(z | x t z))2
) ] (11)
where p(x t z) is given in equation (5) and pmiss(z | x t 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 pmiss.
Appendix F discusses situations where training on incomplete data by EM is possible.

heading: 4 A LOSS FUNCTION AND DECODING METHOD
text: It is sometimes useful to find a single hypothesis ẑ 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:
ẑ = arg minz∈Z ∑ z∗∈Z p(z∗ | x)D(z, z∗) ≈ arg minz∈Z M∑ m=1 wmD(z, zm) (13)
where D(z∗, z) is the distance between z∗ and z. We now propose a specific distance function D.

heading: 4.1 OPTIMAL TRANSPORT DISTANCE
text: 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 Cdelete, and each unaligned event in z∗ incurs an insertion cost of Cinsert. Thus,
D(z∗, z) = min a∈A(z∗,z) D(z∗, z,a) (14)
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., Cinsert = Cdelete = C, so that the distance is symmetric.

heading: 4.2 APPROXIMATE MBR DECODING
text: Since aligned events must have the same type, the MBR problem decomposes into separately choosing a set ẑ(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 ẑ 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:
Theorem 1. Given {zm}Mm=1, if we define zt = ⊔M m=1 zm, then:
∃ẑ ∈ P(zt) such that M∑ m=1 wmD(zm, ẑ) = min z∈Z M∑ m=1 wmD(zm, z) (15)
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 zt 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 t1 in zt. Recall that the alignment cost between these two events is |t− t1|, so we can simply change z∗ by moving t to t1 to minimize its total alignment cost. Then what if it is aligned to events at t1 < t2? While t < t1, increasing t always decreases its total alignment cost |t− t1|+ |t− t2| = t1 + t2−2t. Similarly, we should decrease t if t > t2. When t1 ≤ t ≤ t2, the alignment cost is fixed at t2 − t1. Therefore, the total alignment cost of this event is always minimized if t is moved to t1 or t2. 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:
ẑ = arg minz∈P(zt) M∑ m=1 wmD(zm, z) (16)
which can be approximately solved although it is still NP-hard to exactly solve.
Our algorithm (details in Algorithm 3 of Appendix D) finds ẑ by iteratively (1) finding its optimal alignment am with each zm (called Align Phase) using this method of section 4.1, and then (2) going through the following phases—each of them will update ẑ and decreases the weighted distance∑M m=1 wmD(zm, ẑ,am) which is an upper bound of ∑M m=1 wmD(zm, ẑ):
Move Phase Thanks to Theorem 1, given fixed |ẑ| and {am}Mm=1, ∑M m=1 wmD(zm, ẑ,am) can be
minimized by simply moving each ti in the current ẑ to one of the times, if there is any, which it is aligned to. Otherwise, we keep ti unchanged.
Delete Phase Then we may delete any event in ẑ and its associated element in each am if∑M m=1 wmD(zm, ẑ,am) 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 ẑ. 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 wmD(zm, ẑ,am) by inserting to ẑ some events
that can be aligned, at low cost, to those in zm which would have been left alone otherwise. Thanks again to Theorem 1, such fortunate events can be chosen from ⊔M m=1 zm.
Note that even though we used D(zm, ẑ,am) 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 am. Details of the full algorithm (and its theoretical guarantee) can be found in Appendix D.

heading: 5 EXPERIMENTS
text: 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).

heading: 5.1 DATASETS
text: 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).

heading: 5.2 EVALUATION AND RESULTS
text: 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
3We show the 2 real datasets only. The figures for the 10 synthetic datasets are boringly similar to these.
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).

heading: 6 DISCUSSION
text: 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 t 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 ki@ti pair, not just a discrete type ki—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.
4This dependence actually makes the distance unintuitive. For example, on interval [0, 100), the distance between z1 = 1, 3, 4, 5 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.
5Besides 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).

heading: A SEQUENTIAL MONTE CARLO DETAILS
text: 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 Hm 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 λpk(t) and λ q k(t).
Integral computation Computing the weight wm 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 (ti,j − ti,j−1)λp(t) and (ti,j − ti,j−1)λq(t) at a random t ∼ Unif(ti,j−1, ti,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 λ∗ = fk(maxt g1(t)+ . . .+ maxt gn(t)) where each summand vkdhd(t) = vkd · oid · (2σ(2cd(t))− 1) is upper-bounded by maxc∈{cid,c̄id} vkd · oid · (2σ(2c) − 1). Note that the coefficients wkd may be either positive or negative.
For particle smoothing, we simply have more summands inside fk so λ∗ = fk(maxt g1(t) + . . . + maxt gn(t) + maxt ḡ1(t) + . . . + maxt ḡn̄(t)) where each extra summand ukdh̄d(t) = ukd · oid · (2σ(2cd(t)) − 1) is upper-bounded by maxc∈{cid,c̄id} ukd · oid · (2σ(2c) − 1) and each ukd is the d-th element of vector v>kV (equation (9)). Note that the oid, cid, c̄id of newly added summands ḡ 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 pmiss(Miss = z | Comp = x t 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 w̃m 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 pmiss(Miss = z | Comp = x t z) factors as∏
`∈indices(xtz)
pmiss((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 pmiss 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 rk(t | H(t)) of going missing given the preceding eventsH(t), it is wise to define λqk(t | H(t)) = λ p k(t | H(t)) · rk(t | H(t)).
We include this extra rk factor in our experiments (section 5), although it is not shown in Algorithm 1. It is particularly simple in our experiments, where rk is constant at 1 or 0 depending on
6This 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).
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 pmiss. That is because section 3.2.1 aims to train λqk(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 rk 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 rk. This yields a more expressive and better-factored family of proposal distributions: missingness is now handled by the known rk factor and the BiLSTM does not have to explain it.
Modifying equation (9) in this way is particularly useful in the special case rk = 0 (i.e., event type k is never missing and should not be proposed). There, it enforces the hard constraint that λqk = 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.

heading: B RIGHT-TO-LEFT CONTINUOUS-TIME LSTM
text: 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: h̄(t) = oi (2σ(2c(t))− 1) for t ∈ (ti−1, ti] (18) where the interval (ti−1, ti) has consecutive observations ki−1@ti−1 and ki@ti as endpoints. At ti, the continuous-time LSTM reads ki@ti and updates the current (decayed) hidden cells c(t) to new initial values ci−1, based on the current (decayed) hidden state h̄(ti), as follows:7
ii−1 ← σ (Wiki + Uih(ti) + di) (19a) f i−1 ← σ (Wfki + Ufh(ti) + df) (19b) zi−1 ← 2σ (Wzki + Uzh(ti) + dz)− 1 (19c) oi−1 ← σ (Woki + Uoh(ti) + do) (19d)
ci−1 ← f i−1 c(ti) + ii−1 zi−1 (20a) c̄i−1 ← f̄ i−1 c̄i + īi−1 zi−1 (20b) δi−1 ← f (Wdki + Udh(ti) + dd) (20c)
The vector ki ∈ {0, 1}K is the ith input: a one-hot encoding of the new event ki, with non-zero value only at the entry indexed by ki. Then, c(t) is given by (21), which continues to control h(t) except that i has now decreased by 1).
c(t) def = c̄i−1 + (ci−1 − c̄i−1) exp (−δi−1 (ti − t)) for t ∈ (ti−1, ti] (21)
On the interval [ti−1, ti), c(t) follows an exponential curve that begins at ci−1 (in the sense that limt→t−i c(t) = ci−1) and decays, as time t decreases, toward c̄i−1.

heading: C OPTIMAL TRANSPORT DISTANCE DETAILS
text: 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.
7The 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 ī in equation (20b) are defined analogously to f and i but with different weights.
Algorithm 1 Sequential Monte Carlo — Neural Hawkes Particle Filtering/Smoothing Input: observed seq. x = k0@t0, . . . , kI+1@tI+1 with k0 = BOS, t0 = 0, kI+1 = EOS, tI+1 = T ;
model p; missingness mechanism pmiss; proposal distribution q; number of particles M ; boolean flags smooth and resample
Output: collection {(z1, w1), . . . , (zM , wM )} of weighted particles 1: procedure SEQUENTIALMONTECARLO(x, p, pmiss, q,M, smooth, resample) 2: for m = 1 to M : . init weighted particles (zm, wm). HistoryHm combines zm with a prefix of x 3: zm ← empty sequence; wm ← 1;Hm ← empty stack 4: if smooth : . use particle smoothing 5: F ← empty stack 6: for i = I downto 0 : . stack of all future observed events 7: push ki+1@ti+1 onto F . as we reach these events, we’ll pop from F and push ontoHm (∀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 : 12: DRAWSEGMENT(i,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 {(zm, wm/ ∑M m=1 wm)}Mm=1 .M particles with weights normalized as in equation (7a) 16: procedure RESAMPLE . has access to global variables 17: for m = 1 to M : . often draws multiple copies of good (high-weight) particles, 0 copies of bad ones 18: z̃m ∼ Categorical({zm 7→ wm/ ∑M m=1 wm}Mm=1) 19: for m = 1 to M : 20: zm ← z̃m; wm ← 1 . update particles and their weights 21: procedure DRAWSEGMENT(i,m) . has access to global variables 22: . p gives info to define function λpk(t) def = λk(t | Hm) 23: . q gives info to define function λqk(t) def = λk(t | Hm,F), or λqk(t) = λ p k(t) if F = ignored 24: . these functions consult state of left-to-right LSTM that’s readHm & 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; ki,j ← ki; ti,j ← ti . ready to observe ith 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: wm ← wm · λpki,j (ti,j) exp (− ∫ ti,j t′=ti,j−1
λp(t′)dt′) . new factor in numerator p of (7b) 29: push ki,j@ti,j ontoHm; t← ti,j . event just generated by p now becomes part of history 30: wm ← wm · pmiss((ki,j@ti,j ∈ Miss) = (j > 0) | Hm) . 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, ti+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 ≥ ti+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; ti,j ← t; ki,j ∼ Categorical({k 7→ λqk(t) λq(t)} K k=1); append ki,j@ti,j to zm 40: wm ← wm/ ( λqki,j (ti,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: D(z1, z2) = 0⇔ z1 = z2, (22)
where z1 and z2 are two sequences. If z1 is not identical to z2, 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 z1 and z2 is zero, which means
Algorithm 2 A Dynamic Programming Algorithm to Find Optimal Transport Distance Input: proposal ẑ; reference z∗ Output: optimal transport distance d; alignment a
1: procedure OTD(z∗, ẑ) 2: d← 0; a← empty collection {} 3: for k ← 1 to K : 4: d(k),a(k) ← DYNAMICPROGRAMMING(ẑ(k), z∗(k)) 5: d← d+ d(k); a← a ∪ a(k) 6: return d,a 7: procedure DYNAMICPROGRAMMING(z∗(k), ẑ(k)) 8: Î ← |ẑ(k)|; I∗ ← |z∗(k)| . ẑ(k) = t̂1, . . . , t̂Î 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 î← 1 to Î : . transport reference of length 0 to proposal of length î 12: Dî,0 ← Dî−1,0 + Cinsert . insert t ∗ î = t̂î to reference (and their prefixes are matched) 13: for i∗ ← 1 to I∗ : . transport preference of length i∗ to proposal of length 0 14: D0,j∗ ← D0,j∗ + Cdelete . delete t∗i∗ (and prefixes are matched) 15: for î← 1 to Î : . proposal prefix of length î 16: for i∗ ← 1 to I∗ : . to match reference of length i∗ 17: Dinsert ← Dî−1,i∗ + Cinsert . if an event token at t̂î is inserted to z ∗(k) 18: Ddelete ← Dî,i∗−1 + Cdelete . if the event token at t ∗ i∗ is deleted from z ∗(k) 19: Dmove ← Dî−1,i∗−1 + |t̂̂i − t∗i∗ | . if the event at t ∗ i∗ of z ∗(k) is aligned to event at t̂î of ẑ (k) 20: Dî,i∗ ← min{Dinsert, Ddelete, Dmove} . choose the edit that yields the shortest distance 21: Pî,i∗ ← arg mine∈{insert,delete,move}De . e represents a kind of edition 22: î← Î; i∗ ← I∗;a← empty collection{} 23: while î > 0 and i∗ > 0 : . back trace 24: if Pî,i∗ = insert : . token t ∗ i∗ is deleted.

heading: 25: î← î− 1
text: 26: if Pî,i∗ = delete : . a token at t̂î is inserted 27: i∗ ← i∗ − 1 28: if Pî,i∗ = move : . token t ∗ i∗ is aligned to t̂î

heading: 29: î← î− 1; i∗ ← i∗ − 1

heading: 30: a← a ∪ {(t̂̂i, t∗i∗)}
text: 31: return DÎ,I∗ ,a
they are already matched without any operations, z1 and z2 must be identical, thus the left direction of equation (22) holds.
OTD is symmetric, that is, D(z1, z2) = D(z2, z1), if we set Cinsert = Cdelete. Suppose that a is an alignment between z1 and z2. It’s easy to see that the only difference between D(z1, z2,a) and D(z2, z1,a)
8 is that the insertion and deletion operations are exchanged. For example, if we delete a token ti ∈ z1 when calculating D(z1, z2,a), we should insert a token at ti to z2 when calculating D(z2, z1,a). If we set Cinsert = Cdelete, we have D(z1, z2,a) = D(z2, z1,a), ∀a ∈ A(z1, z2). (23) Therefore, we could obtain
D(z1, z2) = min a∗∈A(z1,z2)
D(z1, z2,a ∗) = min
a∗∈A(z1,z2) D(z2, z1,a
∗) = D(z2, z1). (24)
Finally let’s prove that OTD satisfies triangle inequality, that is: D(z1, z2) +D(z2, z3) ≥ D(z1, z3), (25)
where z1, z2 and z3 are three sequences. This property could be proved intuitively. Suppose that the operations on z1 with minimal costs to make z1 matched to z2 are denoted by o1, o2, . . . , on1 , and
8We abuse the notation a, which we think could represent both the movement from z1 to z2 and from z2 to z1.
those on z2 to make z2 matched to z3 are denoted by o′1, o ′ 2, . . . , o ′ n2 . oi could be a deletion, insertion or movement on a token. To make z1 matched to z3, one possible way, which is not necessarily the optimal, is to do o1, o2, . . . , on1 , o ′ 1, o ′ 2, . . . , o ′ n2 on z1. Since the total cost is the accumulation of the cost of each operation, and the operations on z1 above to make z1 matched to z3 might not be optimal, the triangle inequality equation (25) holds.

heading: D APPROXIMATE MBR DETAILS
text: 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: Theorem 1. Given {zm}Mm=1, if we define zt = ⊔M m=1 zm, then:
∃ẑ ∈ P(zt) such that M∑ m=1 wmD(zm, ẑ) = min z∈Z M∑ m=1 wmD(zm, z) (26)
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 zt 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, M∑ m=1 wmD(zm, z ∗) = min z∈Z M∑ m=1 wmD(zm, z).
If z∗ ∈ P(zt), the proof is done. If not, suppose there exists a token at ti /∈ zt. We use tl ∈ zt (tr ∈ zt) to denote the token in zt that is left (right) and nearest to ti.9 We will show that if we move ti around, as long as ti ∈ [tl, tr], the weighted optimal transport distance, i.e. ∑M m=1 wmD(zm, z
∗), will neither increase nor decrease.
Suppose a∗m = arg minam∈A(zm,z∗) ∑M m=1 wmD(zm, z
∗,am). Let’s use r(t) to indicate the weighted transport distance of z∗ with fixed alignment if we move ti to t, that is,
r(t) def = M∑ m=1 wmD(zm, z ∗(t),a∗m),
where z∗(t) is the sequence z∗ with ti moved to t. Because z∗(ti) is an optimal decode, and a∗m is the optimal alignment for z∗(ti), we should have
r(ti) = 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 ∈ [tl, tr], which guarantees that it will not cross any other tokens in zt.
Since r(t) is linear to t ∈ [rl, tr] and r(t) gets minimized at ti ∈ (tl, tr), we could conclude that r(t) = r(ti) = Const,∀t ∈ [tl, tr].
9We assume that ti is in between two tokens in zt, and our proof could be easily extended to the case for which this condition is not satisfied.
Since r(t) is the upper bound of the weighted optimal transport distance
∑M
m=1 wmD(zm, z
∗(t)), which also gets the same minimal value at ti ∈ (tl, tr) as r(t), we could conclude that
M∑ m=1 wmD(zm, z ∗(t)) = M∑ m=1 wmD(zm, z ∗(ti)) = Const,∀t ∈ [tl, tr].
Therefore we could move token ti to either tl or tr without increasing the Bayes risk. We could do this movement for each ti /∈ zt to get a new decode ẑ ∈ P(zt), which is also an optimal decode.
D.2 ALGORITHM DETAILS
The detailed algorithm is presented in Algorithm 3.

heading: E EXPERIMENTAL DETAILS
text: 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 pmiss(Miss = z | Comp = x t z) = ∏
ki@ti∈z
ρki ∏
ki@ti∈x
(1− ρki), (27)
meaning that each event in the complete stream xt 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 pmiss(Miss = z | Comp = xt z) in equation (3) is not constant. The factor is 1 if z consists of precisely the events in xt z that ought to go missing, and 0 otherwise. In other words, pmiss here ensures only that z should consist only of events with ρk > 0. Recall from Appendix A that pmiss 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
{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. 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 = tI , 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.
10We rebuilt the system in Python following the original Fortran code of Crites & Barto (1996).
E.5 NEW YORK CITY TAXI DATASET DETAILS
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 t0, 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.

heading: F MONTE CARLO EM
text: 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: p(Obs = x) = ∑ z p(Miss = z,Obs = x) = ∑ z p(x t z)pmiss(z | x) (28) 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:
p(x) = ∑ z p(x t z)pmiss(z | x) q(z | x) q(z | x) dz = Ez∼q(z|x)[ p(x t z)pmiss(z | x) q(z | x) ] (29)
Given a finite number M of proposed particles {zm}Mm=1, this expectation can be estimated with empirical average:
p(x) = 1
M M∑ m=1 p(x t zm)pmiss(zm | x) q(zm | x)
(30)
and it is obvious that
log p(x) ≥ 1 M M∑ m=1 (log p(x t zm) + log pmiss(zm | x)− log q(zm | x)) (31)
11It is “almost impossible” to determine from the data whether the MAR assumption holds (Mohan & Pearl, 2018).
where the right-hand-side (RHS) term is the Evidence Lower Bound (ELBO) that we would maximize in order to maximize the log-likelihood.
The MCEM algorithm is composed of two steps:
E(xpectation)-step We train the proposal distribution q(z | x) using the method in section 3.2.1 and then sample M weighted particles from q(z | x) by calling Algorithm 1.
M(aximization)-step We train the neural Hawkes process p(xtz) by maximizing the ELBO (equation (31)).
Note that in the MAR case, pmiss(z | x) is constant of z so the it can be omitted from the formulation (and thus the algorithms). Also note that, for particle filtering, the proposal distribution q(z | x) is only part of p(x t z) so we ought only to propose particles in the E-step.

heading: G MORE EXPERIMENTS ON MNAR DATA
text: 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 pmiss(Miss = z | Comp = x t 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, pmiss 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.
Algorithm 3 Find an approximately minimum-Bayes-risk sequence of events
Input: collection of weighted particles ZM = {(zm, wm)}Mm=1 Output: decode ẑ
1: procedure APPROXMBR(ZM ) 2: ẑ← empty sequence 3: for k = 1 to K : 4: ẑ(k) ← DECODEK({(z(k)m , wm)}Mm=1) . decode for type-k by calling DECODEK 5: ẑ← ẑ t ẑ(k) 6: return ẑ 7: procedure DECODEK(ZM ) 8: . ZM actually means Z(k)M = {(z (k) m , wm)}Mm=1 throughout the procedure; zm is constant
9: z← arg maxz∈{zm}Mm=1 wm . init decode as highest weighted particle and it is global 10: repeat 11: for m = 1 to M : . Align Phase 12: dm,am ← DYNAMICPROGRAM(zm, z) . call method in Algorithm 2; dm,am are global 13: rmin ← ∑ m wmdm . track the risk of current z 14: z, {dm,am}Mm=1 ← MOVE(z, {zm, dm,am}Mm=1) 15: z, {dm,am}Mm=1 ← DELETE(z, {zm, dm,am}Mm=1) 16: z, {dm,am}Mm=1 ← INSERT(z, {zm, dm,am}Mm=1) 17: until ∑M m=1 wmdm = rmin . risk stops decreasing 18: return z 19: procedure MOVE(z, {zm, dm,am}Mm=1) . Move Phase 20: for t in z : 21: for t′ ∈ {t′ : (t′, t) ∈ ⋃M m=1 am} : . may replace t with t′ which is aligned to t 22: (∀m)d′m ← dm 23: for (t′′,m) ∈ {(t′′,m) : (t′′, t) ∈ am,m ∈ {1, . . . ,M}} : 24: d′m ← d′m − |t′′ − t|+ |t′′ − t′| 25: if ∑ m wmd ′ m < ∑ m wmdm : 26: (∀m)dm ← d′m; t← t′ . t move to t′ for lower risk 27: return z, {dm,am}Mm=1 28: procedure DELETE(z, {zm, dm,am}Mm=1) . Delete Phase 29: for t in z : . may delete this event 30: for m = 1 to M : . update each dm 31: if ∃t′ ∈ zm and (t′, t) ∈ am : . find the only, if any, t′ ∈ zm that is aligned to t 32: d′m ← dm + Cdelete − |t′ − t| . dm changes if we delete t and alignment (t′, t) 33: else . otherwise, one insertion must have been made to match t 34: d′m ← dm − Cinsert . so the insertion cost disappears if we delete t 35: if ∑ m wmd ′ m < ∑ m wmdm : 36: delete t from z; (∀m) delete (t′, t) from am; dm ← d′m 37: return z, {dm,am}Mm=1 38: procedure INSERT(z, {zm, dm,am}Mm=1) . Insert Phase 39: for t ∈ {t : t ∈ ⋃M m=1 zm, t /∈ z} : . may insert t′ if it is not in z yet 40: (∀m)a′m ← am 41: for m = 1 to M : . find t′ in zm that is not aligned (thus in z′) and closest to t (arg min) 42: z′ ← {t′ : ∀t, (t′, t) /∈ am} . z′ may be empty, i.e. all in zm are aligned 43: if z′ is not empty and |(t′ ← arg mint′∈zm∩z′{|t
′ − t|})− t| < Cinsert + Cdelete : 44: d′m ← dm − Cdelete + |t′ − t|; add (t′, t) to a′m 45: else 46: d′m ← dm + Cinsert 47: if ∑ m wmd ′ m < ∑ m wmdm : 48: insert t to z; (∀m)dm ← d′m; am ← a′m 49: return z, {dm,am}Mm=1

Review summaries:
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.