Poisson--Gamma Dynamical Systems

Introduction

Sequentially observed count vectors $\boldsymbol{y}^{(1)},\ldots,\boldsymbol{y}^{(T)}$ are the main object of study in many real-world applications, including text analysis, social network analysis, and recommender systems. Count data pose unique statistical and computational challenges when they are high-dimensional, sparse, and overdispersed, as is often the case in real-world applications. For example, when tracking counts of user interactions in a social network, only a tiny fraction of possible edges are ever active, exhibiting bursty periods of activity when they are. Models of such data should exploit this sparsity in order to scale to high dimensions and be robust to overdispersed temporal patterns. In addition to these characteristics, sequentially observed multivariate count data often exhibit complex dependencies within and across time steps. For example, scientific papers about one topic may encourage researchers to write papers about another related topic in the following year. Models of such data should therefore capture the topic structure of individual documents as well as the excitatory relationships between topics.

Many previous approaches to modeling sequentially observed count data rely on the generalized linear modeling framework mccullagh1989generalized to link the observations to a latent Gaussian space—e.g., via the Poisson–lognormal link bulmer1974fitting . Researchers have used this construction to factorize sequentially observed count matrices under a Poisson likelihood, while modeling the temporal structure using well-studied Gaussian techniques blei2006dynamic ; charlin2015dynamic . Most of these previous approaches assume a simple Gaussian state-space model—i.e., $\boldsymbol{\theta}^{(t)}\sim\mathcal{N}(\boldsymbol{\theta}^{(t-1)},\Delta)$—that lacks the expressive transition structure of the LDS; one notable exception is the Poisson linear dynamical system macke2011empirical . In practice, these approaches exhibit prohibitive computational complexity in high dimensions, and the Gaussian assumption may fail to accommodate the burstiness often inherent to real-world count data kleinberg2003bursty .

We present the Poisson–gamma dynamical system (PGDS)—a new dynamical system, based on the gamma–Poisson construction, that supports the expressive transition structure of the LDS. This model naturally handles overdispersed data. We introduce a new Bayesian nonparametric prior to automatically infer the model’s rank. We develop an elegant and efficient algorithm for inferring the parameters of the transition structure that advances recent work on augmentation schemes for inference in negative binomial models zhou12augment-and-conquer and scales with the number of non-zero counts, thus exploiting the sparsity inherent to real-world count data. We examine the way in which the dynamical gamma–Poisson construction propagates information and derive the model’s steady state, which involves the Lambert W function corless1996lambertw . Finally, we use the PGDS to analyze a diverse range of real-world data sets, showing that it exhibits excellent predictive performance on smoothing and forecasting tasks and infers interpretable latent structure, an example of which is depicted in figure 1.

Poisson–Gamma Dynamical Systems

We can represent a data set of $V$-dimensional sequentially observed count vectors $\boldsymbol{y}^{(1)},\ldots,\boldsymbol{y}^{(T)}$ as a $V\times T$ count matrix $Y$. The PGDS models a single count $y_{v}^{(t)}\in\{0,1,\ldots\}$ in this matrix as follows:

where the latent factors $\phi_{vk}$ and $\theta_{k}^{(t)}$ are both positive, and represent the strength of feature $v$ in component $k$ and the strength of component $k$ at time step $t$, respectively. The scaling factor $\delta^{(t)}$ captures the scale of the counts at time step $t$, and therefore obviates the need to rescale the data as a preprocessing step. We refer to the PGDS as stationary if $\delta^{(t)}\!=\!\delta$ for $t=1,\ldots,T$. We can view the feature factors as a $V\times K$ matrix $\Phi$ and the time-step factors as a $T\times K$ matrix $\Theta$. Because we can also collectively view the scaling factors and time-step factors as a $T\times K$ matrix $\Psi$, where element $\psi_{tk}=\delta^{(t)}\,\theta_{k}^{(t)}$, the PGDS is a form of Poisson matrix factorization: $Y\sim\textrm{Pois}(\Phi\,\Psi^{T})$ (canny2004gap, ; cemgil09bayesian, ; zhou2011beta, ; gopalan15scalable, ).

To model the strength of each component, we introduce $K$ component weights $\boldsymbol{\nu}=(\nu_{1},\dots,\nu_{K})$ and place a shrinkage prior over them. We assume that the time-step factors and transition weights for component $k$ are tied to its component weight $\nu_{k}$. Specifically, we define the following structure:

where $\boldsymbol{\pi}_{k}=(\pi_{1k},\ldots,\pi_{Kk})$ is the $k^{\textrm{th}}$ column of $\Pi$. Because $\sum_{k_{1}=1}^{K}\pi_{k_{1}k}=1$, we can interpret $\pi_{k_{1}k}$ as the probability of transitioning from component $k$ to component $k_{1}$. (We note that interpreting $\Pi$ as a stochastic transition matrix relates the PGDS to the discrete hidden Markov model.) For a fixed value of $\gamma_{0}$, increasing $K$ will encourage many of the component weights to be small. A small value of $\nu_{k}$ will shrink $\theta^{(1)}_{k}$, as well as the transition weights in the $k^{\textrm{th}}$ row of $\Pi$. Small values of the transition weights in the $k^{\textrm{th}}$ row of $\Pi$ therefore prevent component $k$ from being excited by the other components and by itself. Specifically, because the shape parameter for the gamma prior over $\theta^{(t)}_{k}$ involves a linear combination of $\boldsymbol{\theta}^{(t-1)}$ and the transition weights in the $k^{\textrm{th}}$ row of $\Pi$, small transition weights will result in a small shape parameter, shrinking $\theta^{(t)}_{k}$. Thus, the component weights play a critical role in the PGDS by enabling it to automatically turn off any unneeded capacity and avoid overfitting.

Finally, we place Dirichlet priors over the feature factors and draw the other parameters from a non-informative gamma prior: $\boldsymbol{\phi}_{k}=(\phi_{1k},\ldots,\phi_{Vk})\sim\textrm{Dir}(\eta_{0},\ldots,\eta_{0})$ and $\delta^{(t)},\xi,\beta\sim\textrm{Gam}(\epsilon_{0},\epsilon_{0})$. The PGDS therefore has four positive hyperparameters to be set by the user: $\tau_{0}$, $\gamma_{0}$, $\eta_{0}$, and $\epsilon_{0}$.

Bayesian nonparametric interpretation: As $K\rightarrow\infty$, the component weights and their corresponding feature factor vectors constitute a draw $G=\sum_{k=1}^{\infty}\nu_{k}\mathds{1}_{\boldsymbol{\phi}_{k}}$ from a gamma process $\textrm{GamP}\left(G_{0},\,\beta\right)$, where $\beta$ is a scale parameter and $G_{0}$ is a finite and continuous base measure over a complete separable metric space $\Omega$ ferguson73bayesian . Models based on the gamma process have an inherent shrinkage mechanism because the number of atoms with weights greater than $\varepsilon>0$ follows a Poisson distribution with a finite mean—specifically, $\textrm{Pois}(\gamma_{0}\int_{\varepsilon}^{\infty}\textrm{d}\nu\ \nu^{-1}\exp{(-\beta\,\nu)})$, where $\gamma_{0}=G_{0}(\Omega)$ is the total mass under the base measure. This interpretation enables us to view the priors over $\Pi$ and $\Theta$ as novel stochastic processes, which we call the column-normalized relational gamma process and the recurrent gamma process, respectively. We provide the definitions of these processes in the supplementary material.

Non-count observations: The PGDS can also model non-count data by linking the observed vectors to latent counts. A binary observation $b_{v}^{(t)}$ can be linked to a latent Poisson count $y_{v}^{(t)}$ via the Bernoulli–Poisson distribution: $b_{v}^{(t)}=\mathds{1}(y_{v}^{(t)}\geq 1)$ and $y_{v}^{(t)}\sim\textrm{Pois}(\delta^{(t)}\sum_{k=1}^{K}\phi_{vk}\,\theta^{(t)}_{k})$ (zhou15infinite, ). Similarly, a real-valued observation $r_{v}^{(t)}$ can be linked to a latent Poisson count $y_{v}^{(t)}$ via the Poisson randomized gamma distribution zhou2015gamma . Finally, Basbug and Engelhardt (basbug2016hierarchical, ) recently showed that many types of non-count matrices can be linked to a latent count matrix via the compound Poisson distribution (adelson1966compound, ).

MCMC Inference

MCMC inference for the PGDS consists of drawing samples of the model parameters from their joint posterior distribution given an observed count matrix $Y$ and the model hyperparameters $\tau_{0}$, $\gamma_{0}$, $\eta_{0}$, $\epsilon_{0}$. In this section, we present a Gibbs sampling algorithm for drawing these samples. At a high level, our approach is similar to that used to develop Gibbs sampling algorithms for several other related models zhou12augment-and-conquer ; zhou15negative ; acharya2015nonparametric ; zhou15infinite ; however, we extend this approach to handle the unique properties of the PGDS. The main technical challenge is sampling $\Theta$ from its conditional posterior, which does not have a closed form. We address this challenge by introducing a set of auxiliary variables. Under this augmented version of the model, marginalizing over $\Theta$ becomes tractable and its conditional posterior has a closed form. Moreover, by introducing these auxiliary variables and marginalizing over $\Theta$, we obtain an alternative model specification that we can subsequently exploit to obtain closed-form conditional posteriors for $\Pi$, $\boldsymbol{\nu}$, and $\xi$. We marginalize over $\Theta$ by performing a “backward filtering” pass, starting with $\boldsymbol{\theta}^{(T)}$. We repeatedly exploit the following three definitions in order to do this.

Definition 1: If $y_{\boldsymbol{\cdot}}\!=\!\sum_{n=1}^{N}y_{n}$, where $y_{n}\sim\textrm{Pois}(\theta_{n})$ are independent Poisson-distributed random variables, then $(y_{1},\ldots,y_{N})\sim\textrm{Mult}(y_{\boldsymbol{\cdot}},(\frac{\theta_{1}}{\sum_{n=1}^{N}\theta_{n}},\ldots,\frac{\theta_{N}}{\sum_{n=1}^{N}\theta_{n}}))$ and $y_{\boldsymbol{\cdot}}\sim\textrm{Pois}(\sum_{n=1}^{N}\theta_{n})$ kingman72poisson ; Dunson05bayesianlatent .

Definition 2: If $y\sim\textrm{Pois}(c\,\theta)$, where $c$ is a constant, and $\theta\sim\textrm{Gam}(a,b)$, then $y\sim\textrm{NB}(a,\frac{c}{b+c})$ is a negative binomial–distributed random variable. We can equivalently parameterize it as $y\sim\textrm{NB}(a,g(\zeta))$, where $g(z)=1-\exp{(-z)}$ is the Bernoulli–Poisson link zhou15infinite and $\zeta=\ln{(1+\frac{c}{b})}$.

Definition 3: If $y\sim\textrm{NB}(a,g(\zeta))$ and $l\sim\textrm{CRT}(y,a)$ is a Chinese restaurant table–distributed random variable, then $y$ and $l$ are equivalently jointly distributed as $y\sim\textrm{SumLog}(l,g(\zeta))$ and $l\sim\textrm{Pois}(a\,\zeta)$ zhou15negative . The sum logarithmic distribution is further defined as the sum of $l$ independent and identically logarithmic-distributed random variables—i.e., $y=\sum_{i=1}^{l}x_{i}$ and $x_{i}\sim\textrm{Log}(g(\zeta))$.

Marginalizing over $\Theta$: We first note that we can re-express the Poisson likelihood in equation 1 in terms of latent subcounts cemgil09bayesian : $y_{v}^{(t)}=y_{v\boldsymbol{\cdot}}^{(t)}=\sum_{k=1}^{K}y_{vk}^{(t)}$ and $y_{vk}^{(t)}\sim\textrm{Pois}(\delta^{(t)}\,\phi_{vk}\,\theta_{k}^{(t)})$. We then define $y_{\boldsymbol{\cdot}k}^{(t)}=\sum_{v=1}^{V}y_{vk}^{(t)}$. Via definition 1, we obtain $y_{\boldsymbol{\cdot}k}^{(t)}\sim\textrm{Pois}(\delta^{(t)}\,\theta_{k}^{(t)})$ because $\sum_{v=1}^{V}\phi_{vk}=1$.

We start with $\theta_{k}^{(T)}$ because none of the other time-step factors depend on it in their priors. Via definition 2, we can immediately marginalize over $\theta_{k}^{(T)}$ to obtain the following equation:

Next, we marginalize over $\theta_{k}^{(T-1)}$. To do this, we introduce an auxiliary variable: $l_{k}^{(T)}\sim\textrm{CRT}(y_{\boldsymbol{\cdot}k}^{(T)},\tau_{0}\textstyle{\sum}_{k_{2}=1}^{K}\pi_{kk_{2}}\,\theta_{k_{2}}^{(T-1)})$. We can then re-express the joint distribution over $y_{\boldsymbol{\cdot}k}^{(T)}$ and $l_{k}^{(T)}$ as

We are still unable to marginalize over $\theta_{k}^{(T-1)}$ because it appears in a sum in the parameter of the Poisson distribution over $l_{k}^{(T)}$; however, via definition 1, we can re-express this distribution as

We then define $l_{\boldsymbol{\cdot}k}^{(T)}=\sum_{k_{1}=1}^{K}l_{k_{1}k}^{(T)}$. Again via definition 1, we can express the distribution over $l_{\boldsymbol{\cdot}k}^{(T)}$ as $l_{\boldsymbol{\cdot}k}^{(T)}\sim\textrm{Pois}(\zeta^{(T)}\,\tau_{0}\,\theta_{k}^{(T-1)})$. We note that this expression does not depend on the transition weights because $\sum_{k_{1}=1}^{K}\pi_{k_{1}k}=1$. We also note that definition 1 implies that $(l_{1k}^{(T)},\ldots,l_{Kk}^{(T)})\sim\textrm{Mult}(l_{\boldsymbol{\cdot}k}^{(T)},(\pi_{1},\ldots,\pi_{K}))$. Next, we introduce $m_{k}^{(T-1)}=y_{\boldsymbol{\cdot}k}^{(T-1)}+l_{\boldsymbol{\cdot}k}^{(T)}$, which summarizes all of the information about the data at time steps $T-1$ and $T$ via $y_{\boldsymbol{\cdot}k}^{(T-1)}$ and $l_{\boldsymbol{\cdot}k}^{(T)}$, respectively. Because $y_{\boldsymbol{\cdot}k}^{(T-1)}$ and $l_{\boldsymbol{\cdot}k}^{(T)}$ are both Poisson distributed, we can use definition 1 to obtain

Combining this likelihood with the gamma prior in equation 1, we can marginalize over $\theta_{k}^{(T-1)}$:

We then introduce $l_{k}^{(T-1)}\sim\textrm{CRT}(m_{k}^{(T-1)},\tau_{0}\sum_{k_{2}=1}^{K}\pi_{kk_{2}}\,\theta_{k_{2}}^{(T-2)})$ and re-express the joint distribution over $l_{k}^{(T-1)}$ and $m_{k}^{(T-1)}$ as the product of a Poisson and a sum logarithmic distribution, similar to equation 4. This then allows us to marginalize over $\theta_{k}^{(T-2)}$ to obtain a negative binomial distribution. We can repeat the same process all the way back to $t=1$, where marginalizing over $\theta^{(1)}_{k}$ yields $m_{k}^{(1)}\sim\textrm{NB}(\tau_{0}\,\nu_{k},g(\zeta^{(1)}))$. We note that just as $m_{k}^{(t)}$ summarizes all of the information about the data at time steps $t,\ldots,T$, $\zeta^{(t)}=\ln{(1+\frac{\delta^{(t)}}{\tau_{0}}+\zeta^{(t+1)})}$ summarizes all of the information about $\delta^{(t)},\ldots,\delta^{(T)}$.

As we mentioned previously, introducing these auxiliary variables and marginalizing over $\Theta$ also enables us to define an alternative model specification that we can exploit to obtain closed-form conditional posteriors for $\Pi$, $\boldsymbol{\nu}$, and $\xi$. We provide part of its generative process in figure 2. We define $m_{k}^{(T)}=y_{\boldsymbol{\cdot}k}^{(T)}+l_{\boldsymbol{\cdot}k}^{(T+1)}$, where $l_{\boldsymbol{\cdot}k}^{(T+1)}=0$, and $\zeta^{(T+1)}=0$ so that we can present the alternative model specification concisely.

Steady state: We draw particular attention to the backward pass $\zeta^{(t)}=\ln{(1+\frac{\delta^{(t)}}{\tau_{0}}+\zeta^{(t+1)})}$ that propagates information about $\delta^{(t)},\dots,\delta^{(T)}$ as we marginalize over $\Theta$. In the case of the stationary PGDS—i.e., $\delta^{(t)}=\delta$—the backward pass has a fixed point that we define in the following proposition.

Gibbs sampling algorithm: Given $Y$ and the hyperparameters, Gibbs sampling involves resampling each auxiliary variable or model parameter from its conditional posterior. Our algorithm involves a “backward filtering” pass and a “forward sampling” pass, which together form a “backward filtering–forward sampling” algorithm. We use $-\setminus\Theta^{(\geq t)}$ to denote everything excluding $\boldsymbol{\theta}^{(t)},\ldots,\boldsymbol{\theta}^{(T)}$.

Sampling the auxiliary variables: This step is the “backward filtering” pass. For the stationary PGDS in its steady state, we first compute $\zeta^{*}$ and draw $(l_{\boldsymbol{\cdot}k}^{(T+1)}\,|\,-)\sim\textrm{Pois}(\zeta^{\star}\,\tau_{0}\,\theta_{k}^{(T)})$. For the other variants of the model, we set $l_{\boldsymbol{\cdot}k}^{(T+1)}=\zeta^{(T+1)}=0$. Then, working backward from $t=T,\ldots,2$, we draw

After using equations 8 and 9 for all $k=1,\ldots,K$, we then set $l_{\boldsymbol{\cdot}k}^{(t)}=\textstyle{\sum}_{k_{1}=1}^{K}l_{k_{1}k}^{(t)}$. For the non-steady-state variants, we also set $\zeta^{(t)}=\ln{(1+\frac{\delta^{(t)}}{\tau_{0}}+\zeta^{(t+1)})}$; for the steady-state variant, we set $\zeta^{(t)}=\zeta^{*}$.

Sampling $\Theta$: We sample $\Theta$ from its conditional posterior by performing a “forward sampling” pass, starting with $\boldsymbol{\theta}^{(1)}$. Conditioned on the values of $l_{\boldsymbol{\cdot}k}^{(2)},\ldots,l_{\boldsymbol{\cdot}k}^{(T+1)}$ and $\zeta^{(2)},\ldots,\zeta^{(T+1)}$ obtained via the “backward filtering” pass, we sample forward from $t=1,\ldots,T$, using the following equations:

Sampling $\Pi$: The alternative model specification, with $\Theta$ marginalized out, assumes that $(l_{1k}^{(t)},\ldots,l_{Kk}^{(t)})\sim\textrm{Mult}(l_{\boldsymbol{\cdot}k}^{(t)},(\pi_{1k},\ldots,\pi_{Kk}))$. Therefore, via Dirichlet–multinomial conjugacy,

Sampling $\boldsymbol{\nu}$ and $\xi$: We use the alternative model specification to obtain closed-form conditional posteriors for $\nu_{k}$ and $\xi$. First, we marginalize over $\boldsymbol{\pi}_{k}$ to obtain a Dirichlet–multinomial distribution. When augmented with a beta-distributed auxiliary variable, the Dirichlet–multinomial distribution is proportional to the negative binomial distribution zhou2016nonparametric . We draw such an auxiliary variable, which we use, along with negative binomial augmentation schemes, to derive closed-form conditional posteriors for $\nu_{k}$ and $\xi$. We provide these posteriors, along with their derivations, in the supplementary material.

We also provide the conditional posteriors for the remaining model parameters—$\Phi$, $\delta^{(1)},\ldots,\delta^{(T)}$, and $\beta$—which we obtain via Dirichlet–multinomial, gamma–Poisson, and gamma–gamma conjugacy.

Experiments

In this section, we compare the predictive performance of the PGDS to that of the LDS and that of gamma process dynamic Poisson factor analysis (GP-DPFA) acharya2015nonparametric . GP-DPFA models a single count in $Y$ as $y_{v}^{(t)}\sim\textrm{Pois}(\sum_{k=1}^{K}\lambda_{k}\,\phi_{vk}\,\theta_{k}^{(t)})$, where each component’s time-step factors evolve as a simple gamma Markov chain, independently of those belonging to the other components: $\theta_{k}^{(t)}\sim\textrm{Gam}(\theta_{k}^{(t-1)},c^{(t)})$. We consider the stationary variants of all three models.We used the pykalman Python library for the LDS and implemented GP-DPFA ourselves. We used five data sets, and tested each model on two time-series prediction tasks: smoothing—i.e., predicting $y_{v}^{(t)}$ given $y_{v}^{(1)},\ldots,y_{v}^{(t-1)},y_{v}^{(t+1)},\ldots,y_{v}^{(T)}$—and forecasting—i.e., predicting $y_{v}^{(T+s)}$ given $y_{v}^{(1)},\ldots,y_{v}^{(T)}$ for some $s\in\{1,2,\ldots\}$ durbin2012time . We provide brief descriptions of the data sets below before reporting results.

Global Database of Events, Language, and Tone (GDELT): GDELT is an international relations data set consisting of country-to-country interaction events of the form “country $i$ took action $a$ toward country $j$ at time $t$,” extracted from news corpora. We created five count matrices, one for each year from 2001 through 2005. We treated directed pairs of countries $i{\rightarrow}j$ as features and counted the number of events for each pair during each day. We discarded all pairs with fewer than twenty-five total events, leaving $T=365$, around $V\approx 9,000$, and three to six million events for each matrix.

Integrated Crisis Early Warning System (ICEWS): ICEWS is another international relations event data set extracted from news corpora. It is more highly curated than GDELT and contains fewer events. We therefore treated undirected pairs of countries $i{\leftrightarrow}j$ as features. We created three count matrices, one for 2001–2003, one for 2004–2006, and one for 2007–2009. We counted the number of events for each pair during each three-day time step, and again discarded all pairs with fewer than twenty-five total events, leaving $T=365$, around $V\approx 3,000$, and 1.3 to 1.5 million events for each matrix.

State-of-the-Union transcripts (SOTU): The SOTU corpus contains the text of the annual SOTU speech transcripts from 1790 through 2014. We created a single count matrix with one column per year. After discarding stopwords, we were left with $T=225$, $V=7,518$, and 656,949 tokens.

DBLP conference abstracts (DBLP): DBLP is a database of computer science research papers. We used the subset of this corpus that Acharya et al. used to evaluate GP-DPFA acharya2015nonparametric . This subset corresponds to a count matrix with $T=14$ columns, $V=1,771$ unique word types, and 13,431 tokens.

NIPS corpus (NIPS): The NIPS corpus contains the text of every NIPS conference paper from 1987 to 2003. We created a single count matrix with one column per year. We treated unique word types as features and discarded all stopwords, leaving $T=17$, $V=9,836$, and 3.1 million tokens.

Experimental design: For each matrix, we created four masks indicating some randomly selected subset of columns to treat as held-out data. For the event count matrices, we held out six (non-contiguous) time steps between $t=2$ and $t=T-3$ to test the models’ smoothing performance, as well as the last two time steps to test their forecasting performance. The other matrices have fewer time steps. For the SOTU matrix, we therefore held out five time steps between $t=2$ and $t=T-2$, as well as $t=T$. For the NIPS and DBLP matrices, which contain substantially fewer time steps than the SOTU matrix, we held out three time steps between $t=2$ and $t=T-2$, as well as $t=T$.

For each matrix, mask, and model combination, we ran inference four times.For the PGDS and GP-DPFA we used $K=100$. For the PGDS, we set $\tau_{0}=1$, $\gamma_{0}=50$, $\eta_{0}=\epsilon_{0}=0.1$. We set the hyperparameters of GP-DPFA to the values used by Acharya et al. acharya2015nonparametric . For the LDS, we used the default hyperparameters for pykalman, and report results for the best-performing value of $K\in\{5,10,25,50\}$. For the PGDS and GP-DPFA, we performed 6,000 Gibbs sampling iterations, imputing the missing counts from the “smoothing” columns at the same time as sampling the model parameters. We then discarded the first 4,000 samples and retained every hundredth sample thereafter. We used each of these samples to predict the missing counts from the “forecasting” columns. We then averaged the predictions over the samples. For the LDS, we ran EM to learn the model parameters. Then, given these parameter values, we used the Kalman filter and smoother kalman1960new to predict the held-out data. In practice, for all five data sets, $V$ was too large for us to run inference for the LDS, which is $\mathcal{O}((K+V)^{3})$ ghahramani1999learning , using all $V$ features. We therefore report results from two independent sets of experiments: one comparing all three models using only the top $V=1,000$ features for each data set, and one comparing the PGDS to just GP-DPFA using all the features. The first set of experiments is generous to the LDS because the Poisson distribution is well approximated by the Gaussian distribution when its mean is large.

Results: We used two error measures—mean relative error (MRE) and mean absolute error (MAE)—to compute the models’ smoothing and forecasting scores for each matrix and mask combination. We then averaged these scores over the masks. For the data sets with multiple matrices, we also averaged the scores over the matrices. The two error measures differ as follows: MRE accommodates the scale of the data, while MAE does not. This is because relative error—which we define as $\tfrac{|y_{v}^{(t)}-\hat{y}_{v}^{(t)}|}{1+y_{v}^{(t)}}$, where $y_{v}^{(t)}$ is the true count and $\hat{y}_{v}^{(t)}$ is the prediction—divides the absolute error by the true count and thus penalizes overpredictions more harshly than underpredictions. MRE is therefore an especially natural choice for data sets that are bursty—i.e., data sets that exhibit short periods of activity that far exceed their mean. Models that are robust to these kinds of overdispersed temporal patterns are less likely to make overpredictions following a burst, and are therefore rewarded accordingly by MRE.

In table 1, we report the MRE and MAE scores for the experiments using the top $V=1,000$ features. We also report the average burstiness of each data set. We define the burstiness of feature $v$ in matrix $Y$ to be $\hat{B}_{v}=\frac{1}{T-1}\sum_{t=1}^{T-1}\frac{|y_{v}^{(t+1)}-y_{v}^{(t)}|}{\hat{\mu}_{v}}$, where $\hat{\mu}_{v}=\frac{1}{T}\sum_{t=1}^{T}y_{v}^{(t)}$. For each data set, we calculated the burstiness of each feature in each matrix, and then averaged these values to obtain an average burstiness score $\hat{B}$. The PGDS outperformed the LDS and GP-DPFA on seven of the ten prediction tasks when we used MRE to measure the models’ performance; when we used MAE, the PGDS outperformed the other models on five of the tasks. In the supplementary material, we also report the results for the experiments comparing the PGDS to GP-DPFA using all the features. The superiority of the PGDS over GP-DPFA is even more pronounced in these results. We hypothesize that the difference between these models is related to the burstiness of the data. For both error measures, the only data set for which GP-DPFA outperformed the PGDS on both tasks was the NIPS data set. This data set has a substantially lower average burstiness score than the other data sets. We provide visual evidence in figure 3, where we display $y_{v}^{(t)}$ over time for the top four features in the NIPS and ICEWS data sets. For the former, the features evolve smoothly; for the latter, they exhibit bursts of activity.

Exploratory analysis: We also explored the latent structure inferred by the PGDS. Because its parameters are positive, they are easy to interpret. In figure 1, we depict three components inferred from the NIPS data set. By examining the time-step factors and feature factors for these components, we see that they capture the decline of research on neural networks between 1987 and 2003, as well as the rise of Bayesian methods in machine learning. These patterns match our prior knowledge.

In figure 4, we depict the three components with the largest component weights inferred by the PGDS from the 2003 GDELT matrix. The top component is in blue, the second is in green, and the third is in red. For each component, we also list the sixteen features (directed pairs of countries) with the largest feature factors. The top component (blue) is most active in March and April, 2003. Its features involve USA, Iraq (IRQ), Great Britain (GBR), Turkey (TUR), and Iran (IRN), among others. This component corresponds to the 2003 invasion of Iraq. The second component (green) exhibits a noticeable increase in activity immediately after April, 2003. Its top features involve Israel (ISR), Palestine (PSE), USA, and Afghanistan (AFG). The third component exhibits a large burst of activity in August, 2003, but is otherwise inactive. Its top features involve North Korea (PRK), South Korea (KOR), Japan (JPN), China (CHN), Russia (RUS), and USA. This component corresponds to the six-party talks—a series of negotiations between these six countries for the purpose of dismantling North Korea’s nuclear program. The first round of talks occurred during August 27–29, 2003.

In figure 5, we also show the component weights for the top ten components, along with the corresponding subset of the transition matrix $\Pi$. There are two components with weights greater than one: the components that are depicted in blue and green in figure 4. The transition weights in the corresponding rows of $\Pi$ are also large, meaning that other components are likely to transition to them. As we mentioned previously, the GDELT data set was extracted from news corpora. Therefore, patterns in the data primarily reflect patterns in media coverage of international affairs. We therefore interpret the latent structure inferred by the PGDS in the following way: in 2003, the media briefly covered various major events, including the six-party talks, before quickly returning to a backdrop of the ongoing Iraq war and Israeli–Palestinian relations. By inferring the kind of transition structure depicted in figure 5, the PGDS is able to model persistent, long-term temporal patterns while accommodating the burstiness often inherent to real-world count data. This ability is what enables the PGDS to achieve superior predictive performance over the LDS and GP-DPFA.

Summary

We introduced the Poisson–gamma dynamical system (PGDS)—a new Bayesian nonparametric model for sequentially observed multivariate count data. This model supports the expressive transition structure of the linear dynamical system, and naturally handles overdispersed data. We presented a novel MCMC inference algorithm that remains efficient for high-dimensional data sets, advancing recent work on augmentation schemes for inference in negative binomial models. Finally, we used the PGDS to analyze five real-world data sets, demonstrating that it exhibits superior smoothing and forecasting performance over two baseline models and infers highly interpretable latent structure.

We thank David Belanger, Roy Adams, Kostis Gourgoulias, Ben Marlin, Dan Sheldon, and Tim Vieira for many helpful conversations. This work was supported in part by the UMass Amherst CIIR and in part by NSF grants SBE-0965436 and IIS-1320219. Any opinions, findings, conclusions, or recommendations are those of the authors and do not necessarily reflect those of the sponsors.

References