# Online Bayesian changepoint detection for network Poisson processes with
community structure

Joshua Corneck  Department of Mathematics, Imperial College London
180 Queen’s Gate, SW7 2AZ, London  Edward A. K. Cohen  Department of
Mathematics, Imperial College London
180 Queen’s Gate, SW7 2AZ, London  James S. Martin  Department of Mathematics,
Imperial College London
180 Queen’s Gate, SW7 2AZ, London  Francesco Sanna Passino  Department of
Mathematics, Imperial College London
180 Queen’s Gate, SW7 2AZ, London

###### Abstract

Network point processes often exhibit latent structure that govern the
behaviour of the sub-processes. It is not always reasonable to assume that
this latent structure is static, and detecting when and how this driving
structure changes is often of interest. In this paper, we introduce a novel
online methodology for detecting changes within the latent structure of a
network point process. We focus on block-homogeneous Poisson processes, where
latent node memberships determine the rates of the edge processes. We propose
a scalable variational procedure which can be applied on large networks in an
online fashion via a Bayesian forgetting factor applied to sequential
variational approximations to the posterior distribution. The proposed
framework is tested on simulated and real-world data, and it rapidly and
accurately detects changes to the latent edge process rates, and to the latent
node group memberships, both in an online manner. In particular, in an
application on the Santander Cycles bike-sharing network in central London, we
detect changes within the network related to holiday periods and lockdown
restrictions between 2019 and 2020.

Keywords — network point process, online variational inference, stochastic
blockmodel, streaming data.

##  1 Introduction

Network data describe complex relationships between a large number of
entities, called nodes, linked via edges. In particular, much of such data is
continuous-time valued and dynamic in nature, with nodes often exhibiting a
latent community structure that drives interactions (Amini et al., 2013). By
allowing temporal point processes to exist on the edges of the network, and
not solely focusing on the adjacency matrix, such interaction data can be more
accurately understood. Common examples of these network structures include
email exchanges within a company (Klimt and Yang, 2004) or posts on a social
media platform (Paranjape et al., 2017), where employer hierarchy or political
leanings, respectively, could drive connections between individuals.
Furthermore, it is not always reasonable to assume that such structure is
static. In the case of email exchanges, a change in an employee’s role within
their organisation would alter the latent community structure of the network.
Similarly, in social networks, a sporting or political event would likely
affect the frequency of interactions only between users with particular
characteristics and topics of interest. Detecting such changepoints is crucial
in many real-world applications, such as in the context of computer networks,
where modifications of network structure can signify malicious behaviour (see,
for example, Hallgren et al., 2023). The detection of such changes in an
online manner will be the focus of this paper.

Much of the original work on network models focused on simple, binary
connections between nodes. One of the most common of these modelling
approaches is the Stochastic Block Model (SBM; Holland et al., 1983), which
models connections between nodes via a latent community structure on a
network. This model has undergone significant development, with theoretical
results extending to include node degree heterogeneity (Karrer and Newman,
2011; Amini et al., 2013; Gao et al., 2018), valued edges (Mariadassou et al.,
2010), dynamic community structure in discrete time (Matias and Miele, 2017;
Yang et al., 2011; Xu and Hero, 2014), and missing data (Mariadassou and
Tabouy, 2020); refer to Lee and Wilkinson (2019) for a comprehensive review.
In spite of their simplicity, SBMs can be used to consistently estimate and
approximate much more complex network objects, such as exchangeable random
graphs and their associated graphons, and thus provide a good approximation
for any exchangable random graph model (Airoldi et al., 2013). However, most
of these advancements are based on aggregated, discrete time data.

There are several statistical methods for analysing dynamic networks, normally
defined as graphs with time-evolving edge-connections. Interested readers are
referred to Holme (2015) for a review. However, most of this methodology does
not capture instantaneous interactions, but relies instead on aggregations
(see, for example, Shlomovich et al., 2022). To combat this, Matias et al.
(2018) developed a model for recurrent interaction events, which they call the
semiparametric stochastic block model. In this framework, interactions are
modelled by a conditional inhomogeneous Poisson process, with intensities
depending solely on the latent group memberships of the interacting
individuals. Perry and Wolfe (2013) uses a Cox multiplicative intensity model
with covariates to model the point processes observed on each edge. Sanna
Passino and Heard (2023) models the edge-specific processes via mutually
exciting processes with intensities depending only on node-specific
parameters. These methodologies handle complex temporal data, but they are
offline methods, and work only in the case of a static latent network
structure. Attempts have been made at capturing networks with a dynamic latent
structure, but attention has focused on discrete time data, such as binary or
weighted edges. Matias and Miele (2017) proposed a method for frequentist
inference for a model which extends the SBM to allow for dynamic group
memberships. In that work, the dynamics of the groups are modelled by a
discrete time Markov chain. Heard et al. (2010) developed a two-stage offline
method for anomaly detection in dynamic graphs. However, to the best of our
knowledge, there currently exists no methodology for online detection of
changepoints in the context of continuous-time network data.

The recent work of Fang et al. (2024) represents the first attempt at online
estimation and community detection of a network point process. The authors
build upon the foundation laid by Matias et al. (2018) and extend it to an
online setting, but maintain the assumption of a static latent structure.
While Fang et al. (2024) do offer suggestions as to how their framework could
be extended to incorporate latent dynamics, the methodology developed therein
requires knowing both the adjacency matrix and the number of latent groups a
priori.

In this work, we propose a novel online Bayesian changepoint detection
algorithm for network point processes with a latent community structure among
the nodes. Our methodology is based on utilising forgetting factors within a
Bayesian context to sequentially update the variational approximation to the
posterior distribution of the model parameters when new data are observed
within a stream. We focus on what we will refer to as a dynamic Bayesian
block-homogeneous network Poisson process, where dynamic refers to the fact
that the latent structure of the network is time-varying. As an added benefit,
our method is able to accurately infer the community structure and obtain a
piecewise recovery of the conditional intensity function when we consider
inhomogeneous Poisson processes on the edges. Extensions to our method are
also proposed to simultaneously infer both the adjacency matrix of the network
and the number of latent groups in an online manner, and handle cases where
new groups are created or existing groups are merged into one another.

The remainder of this article is structured as follows: Section 2 describes
the dynamic Bayesian block-homogeneous Poisson process model used in this
work, and the possible local and global changes to the network structure
occurring on such a network model. Section 3 discusses the proposed online
variational Bayesian inference approach via Bayesian forgetting factor, which
is used to sequentially approximate the posterior distribution on the stream.
The performance of the proposed inferential procedure is then tested in
Section 4 on simulated data, and on real-world data from the Santander Cycles
bike-sharing network in Section 5, followed by a discussion.

##  2 Bayesian dynamic block homogeneous Poisson process

###  2.1 The model

We consider a stochastic process on a network, which produces dyadic
interaction data between a set of $N$ nodes observed over time. Let
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a graph, where the set
$\mathcal{V}:=\\{1,\dots,N\\}$ corresponds to nodes, whereas the edge set
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ contains pairs representing
interactions between nodes in $\mathcal{V}$. In particular, we write
$(i,j)\in\mathcal{E}$ if there is a connection from node $i\in\mathcal{V}$ to
node $j\in\mathcal{V}$. Furthermore, we denote the set of all possible edges
between nodes as $\mathcal{R}=\mathcal{V}\times\mathcal{V}$. The graph
$\mathcal{G}$ could be equivalently represented via its adjacency matrix
$\bm{A}=\\{a_{ij}\\}\in\mathbb{\\{}0,1\\}^{N\times N}$ such that
$a_{ij}=\mathbb{I}_{\mathcal{E}}\\{(i,j)\\}$, where
$\mathbb{I}_{\cdot}\\{\cdot\\}$ denotes the indicator function.

We assume that each node $i\in\mathcal{V}$ belongs to a group
$z_{i}\in\mathcal{K}$, where $\mathcal{K}:=\\{1,\dots,K\\}$ and we let
$z=(z_{1},\dots,z_{N})\in\mathcal{K}^{N}$ be a priori independent and
identically distributed with $z_{i}\sim\mathrm{Categorical}(\pi)$ where
$\pi=(\pi_{1},\dots,\pi_{K})$, such that $\pi_{k}\geq 0$ for all
$k\in\mathcal{K}$, and $\sum_{k=1}^{K}\pi_{k}=1$. We also associate to $z_{i}$
a binary vector $(\tilde{z}_{i1},\dots,\tilde{z}_{iK})\in\\{0,1\\}^{K}$, with
$\tilde{z}_{ik}=\mathbb{I}_{\\{k\\}}(z_{i})$, and we write
$\tilde{\bm{Z}}=\\{\tilde{z}_{ik}\\}_{i\in\mathcal{V},k\in\mathcal{K}}\in\\{0,1\\}^{N\times
K}$ to denote the matrix of group memberships. Under a Bayesian framework, we
place a Dirichlet prior distribution on $\pi$ with parameter
$\gamma^{0}\in\mathbb{R}^{K}_{+}$.

We observe a marked point process on the network, consisting of a stream of
triplets $(i_{\ell},j_{\ell},t_{\ell})\in\mathcal{E}\times\mathbb{R}_{+},\
\ell=1,2,\dots$, denoting directed interactions from node $i_{\ell}$ to node
$j_{\ell}$ at time $t_{\ell}$, where $t_{\ell}\leq t_{\ell^{\prime}}$ for
$\ell<\ell^{\prime}$. The associated edge-specific counting process is denoted
by $x_{ij}(\cdot)$, where

$x_{ij}((0,t])\equiv
x_{ij}(t):=\sum_{\ell}\mathbb{I}_{\\{(i,j)\\}}\\{(i_{\ell},j_{\ell})\\}\mathbb{I}_{(0,t]}(t_{\ell}),\quad(i,j)\in\mathcal{E}.$
(2.1)

If $(i,j)\notin\mathcal{E}$, $x_{ij}(t)=0$ for the entire observation period.
In this work, we model the counting process $x_{ij}(\cdot)$ as a Poisson
process with rate $\lambda_{z_{i}z_{j}}\in\mathbb{R}_{+}$, conditional on
$(i,j)\in\mathcal{E}$ and on the node group memberships $z_{i}$ and $z_{j}$.
Also, we place independent conjugate gamma prior distributions on the event
rates $\lambda=\\{\lambda_{km}\\}\in\mathbb{R}_{+}^{K\times K}$. In summary,
the full model can be expressed as follows:

$\displaystyle x_{ij}(t)\mid z_{i},z_{j},\lambda_{z_{i}z_{j}}$
$\displaystyle\sim\mathrm{Poisson}(\lambda_{z_{i}z_{j}}t),$
$\displaystyle\text{for all }(i,j)\in\mathcal{E},$ (2.2)
$\displaystyle\lambda_{km}$
$\displaystyle\sim\mathrm{Gamma}\left(\alpha_{km}^{0},\beta_{km}^{0}\right),$
$\displaystyle\text{for all }k,m\in\mathcal{K},$ (2.3) $\displaystyle
z_{i}\mid\pi$ $\displaystyle\sim\mathrm{Categorical}(\pi),$
$\displaystyle\text{for all }i\in\mathcal{V},$ (2.4) $\displaystyle\pi$
$\displaystyle\sim\mathrm{Dirichlet}\left(\gamma^{0}\right),$ (2.5)

where
$\gamma^{0}=\left(\gamma_{1}^{0},\dots,\gamma_{K}^{0}\right)\in\mathbb{R}^{K}_{+}$,
$\alpha^{0}=\left\\{\alpha^{0}_{km}\right\\}_{k,m\in\mathcal{K}}\in\mathbb{R}^{K\times
K}_{+}$ and
$\beta^{0}=\left\\{\beta_{km}^{0}\right\\}_{k,m\in\mathcal{K}}\in\mathbb{R}^{K\times
K}_{+}$. Note that in (2.2), we use that for a homogeneous Poisson process
with rate $\lambda^{\ast}$, the number of counts for that process on an
interval of length $t$ is Poisson distributed with rate $\lambda^{\ast}t$. The
model in Equation (2.2) is known as block homogeneous Poisson process (BHPP;
Fang et al., 2024), and it is represented in graphical model form in Figure 1.

Figure 1: A directed acyclic graphical model representation of the model in
Equations (2.2)-(2.5).

Under the BHPP, the expected waiting time between events across the entire
network is $\mathcal{O}(N^{-1})$, as described in the following proposition.

######  Proposition 1.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a directed graph with no
isolated nodes. If on each edge lives a homogeneous Poisson process with rate
$\lambda_{z_{i}z_{j}}\in\mathbb{R}_{+}$, only dependent upon node memberships
$z_{i},z_{j}\in\mathcal{K}$, then the expected waiting time between arrivals
for the full network counting process
$x(t)=\sum_{(i,j)\in\mathcal{E}}x_{ij}(t)$ is
$\mathcal{O}\left(N^{-1}\right)$.

###### Proof.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with $|\mathcal{V}|=N$. Suppose
that $\mathcal{G}$ has $M$ connected components. As we have no isolated nodes,
the number of edges is minimised with $M=N/2$ and $|\mathcal{E}|=N/2$ if $N$
is even, or $M=(N-1)/2$ and $|\mathcal{E}|=(N+1)/2$, if $N$ is odd. Without
loss of generality, assume
$\lambda_{11}=\min_{k,m\in\mathcal{K}}\\{\lambda_{km}\\}$ and call the rate of
the full network Poisson process $\lambda$. By superposition, it follows
$\lambda\geq N\lambda_{11}/2$ or $\lambda\geq(N+1)\lambda_{11}/2$, for $n$
even and odd, respectively, and so the expected waiting time $w$, is bounded
as $w\leq 2/\lambda_{11}N$. ∎

The aim of the present paper is to provide an online algorithm for detecting
changepoints; such an algorithm must be able to perform inference at each
iteration before the arrival of any new data. A consequence of Proposition 1
is that for networks of increasing size, this becomes infeasible for large $N$
when observing all events as they arrive, as an online algorithm would need to
be at least as fast. This motivates a further assumption that data arrives as
a stream of batched counts rather than a stream of continuous-valued event
times. Let the batches be observed at times $L_{r},\ r=1,2,\dots$, where
$L_{r}<L_{r^{\prime}}$ for $r<r^{\prime}$. Each data batch takes the form

$\mathcal{D}_{r}:=\left\\{\left(i_{\ell},j_{\ell},t_{\ell}\right):L_{r-1}<t_{\ell}\leq
L_{r}\right\\};\quad r=1,2,\dots;\quad L_{0}\equiv 0.$ (2.6)

We consider the case of a constant time $\Delta$ between batches, so that
$L_{r}-L_{r-1}=\Delta$ for all $r$. Also, we denote the complete set of data
available at time $L_{r}$ as
$\mathcal{D}_{1:r}:=\\{\mathcal{D}_{1},\dots,\mathcal{D}_{r}\\}$. This setup
is illustrated in Figure 2.

Figure 2: Illustration of the setup in Section 2.1 for a network with $N=3$
and $\mathcal{E}=\\{(1,3),(2,1),(2,3)\\}$. Circles represent event times,
coloured according to their edge. Arrival times to each edge are collated into
$\mathcal{D}_{r}$.

Initially, we assume that both the adjacency matrix $\bm{A}$ and the number of
latent groups $K$ are known a priori, similar to the setup in Fang et al.
(2024). However, we note that assuming this knowledge is limiting in practical
applications, since these quantities are usually unknown. Therefore, Sections
3.5 and 3.6 discuss extensions of the BHPP model and corresponding inference
procedures with an unknown adjacency matrix and an unknown number of groups.

###  2.2 Changes to the network

In real-world applications, it is unrealistic to assume that the latent group
memberships $z=(z_{1},\dots,z_{N})$ and rate matrix
$\lambda=\\{\lambda_{km}\\}_{k,m\in\mathcal{K}}$ are constant across the
entire observation period. Our objective is therefore to develop a framework
which could be used to detect global and local changes to the model
parameters.

For local changes, we consider modifications to the latent group structure,
allowing for any $n\leq N$ nodes to change their group membership at any time
$t^{\prime}$. In contrast to the approach of Matias and Miele (2017), in which
the membership of node $i$ over the observation window behaves according to a
discrete-time, irreducible and aperiodic Markov chain, we place no model on
how or when nodes change and allow for any (finite) number of changes in a
given interval. Furthermore, we allow for changes in continuous time, again in
contrast to the discrete-time formulation of Matias and Miele (2017). Our main
objective is only to detect if and when nodes have changed memberships, not to
estimate the underlying group membership switching process. Figure 3
illustrates the effect of a single node switch in the case of $N=4$ nodes.

(a) Observed graph structure

(b) Latent group structure

(c) Change in latent group structure

Figure 3:  Visualisation of an observed graph $\mathcal{G}$ with $N=4$, with a
change in the latent group structure. The changing node is coloured in grey,
and edges are coloured according to the groups they connect.

For global changes, we consider jumps in the intensity between any groups
$k,m\in\mathcal{K}$ at some time $t^{\prime}\in\mathbb{R}_{+}$ of the form
$\lambda_{km}\mapsto\lambda_{km}^{\prime}$
($\lambda_{km}\neq\lambda_{km}^{\prime}$). Again, we are interested in
assessing if the block-specific intensity $\lambda_{km}$ has changed, and we
do not estimate or posit assumptions on the mechanism that leads to such
changes.

##  3 An online variational Bayesian estimation procedure

We present an online inference procedure for tracking the time-evolving latent
structure and parameters of the BHPP model described in Section 2.1. We focus
first on the set-up where the full edge set $\mathcal{E}$ and the number of
latent groups $K$ is known, as in Fang et al. (2024). Next, we consider when
$\mathcal{E}$ is unknown, and then we move to the case of an unknown number of
groups, using a Bayesian nonparametric approach.

###  3.1 Variational Bayesian approximation

Variational Bayesian (VB) inference has the objective to approximate a
posterior distribution when it is not analytically tractable (Wang and Blei,
2019), offering a faster estimation approach when compared to Markov Chain
Monte Carlo methods (Blei et al., 2017). This makes it more suited to the
requirements of an online learning framework. We adopt the terminology of Blei
et al. (2017) in distinguishing between local latent variables
$z_{1},\dots,z_{N}$, that scale with the number of nodes, and global latent
variables $\theta=(\pi,\lambda)$, whose dimension, written $S\in\mathbb{N}$,
does not change with $N$.

In a variational approach, one posits a family of distributions
$\mathcal{F}^{N+S}$ on the parameter space $(\theta,z)$, and seeks to select
the component $q^{\ast}(\theta,z)\in\mathcal{F}^{N+S}$ closest to the true
posterior $p(\theta,z\mid x)$ in the sense of the Kullback-Leibler (KL)
divergence, where $x$ denotes the observed data. A mathematically convenient
choice for $\mathcal{F}^{N+S}$ is the mean-field variational family, the set
of factorisable distributions of the form

$\mathcal{F}^{N+S}=\left\\{q:q(\theta,z)=\prod_{i=1}^{N}q_{z_{i}}(z_{i})\prod_{s=1}^{S}q_{\theta_{s}}(\theta_{s})\right\\}.$
(3.1)

The approximating distribution $q^{*}(\theta,z)\in\mathcal{F}^{N+S}$ is then
given as

$q^{\ast}(\theta,z)=\operatorname*{argmin}_{q\in\mathcal{F}^{N+S}}\mathrm{KL}\left[q(\theta,z)\
\|\ p(\theta,z|x)\right].$ (3.2)

Despite every $q(\theta,z)\in\mathcal{F}^{N+S}$ taking a simple product
structure, the solution $q^{\ast}(\theta,z)$ to (3.2) is usually not available
analytically. A popular method for approximating the global minimum
$q^{\ast}(\theta,z)$ of the KL-divergence is the Coordinate Ascent Variational
Inference (CAVI) algorithm (Bishop, 2006; Blei et al., 2017), which
iteratively updates each component of $q(\theta,z)$ while keeping the others
fixed. Letting $\phi$ denote the $(N+S)$-dimensional vector of latent
variables $\phi=(z,\theta)$, CAVI proceeds as follows:

$\hat{q}_{\phi_{j}}^{(t+1)}(\phi_{j})=\operatorname*{argmin}_{q_{\phi_{j}}}\mathrm{KL}\left[q_{\phi_{j}}(\phi_{j})\prod_{i<j}\hat{q}_{\phi_{i}}^{(t+1)}(\phi_{i})\prod_{i>j}\hat{q}_{\phi_{i}}^{(t)}(\phi_{i})\
\Bigg{|}\Bigg{|}\ p(\phi|x)\right],$ (3.3)

for $j\in\\{1,\dots,N+S\\},\ t=0,1,2,\dots$ until some convergence criterion
is met. Under the mean-field variational approximation, (3.3) has the
following solution (Bishop, 2006):

$\hat{q}_{\phi_{j}}\propto\exp\\{\mathbb{E}_{-\phi_{j}}\left[\log
p(x,\phi)\right]\\},$ (3.4)

where the notation $\mathbb{E}_{-\phi_{j}}$ used to denote the expectation
with respect to all components of $\phi$ except $\phi_{j}$ (Blei et al.,
2017). As per convention with variational inference algorithms, we will drop
the subscript on $q_{\phi_{j}}(\phi_{j})$ for notational convenience, and take
the argument of the variational component to index it (see, for example, Blei
et al., 2017).

It should be noted that CAVI is only guaranteed to achieve a local minimum,
and hence is sensitive to intialisation (Zhang and Zhou, 2020; Blei et al.,
2017). Furthermore, with $\phi=(z,\theta)$, updates in (3.4) have an explicit
ordering. The ordering of the updates can also have implications on the
convergence properties of the algorithm (Ray and Szabó, 2022), similarly to
Gibbs sampling (see, for example, van Dyk and Park, 2008).

###  3.2 Online VB for the dynamic BHPP

Upon receiving the latest batch $\mathcal{D}_{r}$ at time $L_{r}$, the aim of
the online algorithm is to update the estimates of the BHPP parameters based
on the entire history $\mathcal{D}_{1:r}$ of the process. To ensure finite and
constant complexity and memory, a truly online algorithm should be single pass
(see, for example Bifet et al., 2018): each data batch should be inspected
only once, and summaries of previous data batches $\mathcal{D}_{1:(r-1)}$ are
used in conjunction with $\mathcal{D}_{r}$ to update parameter estimates.

The prior-posterior distribution updates within a Bayesian framework offers a
natural way to propagate information forward in time. The posterior
distribution at time $L_{r-1}$ contains the information that was learned from
$\mathcal{D}_{1:(r-1)}$. A naive approach therefore is to pass this through as
the new prior distribution for the update at time $L_{r}$. However, this will
perform poorly in the case of changing latent structure. As the number of
batches increases, the new prior distributions become more concentrated around
the current best estimates. This results in an updated posterior distribution
that is largely dominated by its prior distribution (corresponding to the
posterior distribution up to the previous batch) and which is increasingly
insensitive to new data and changes in the data generating process.

We therefore propose to flatten the posterior obtained at $L_{r-1}$, via
temperature parameters, at the point of passing it through as the prior
distribution for the update at $L_{r}$. This flattening step is in effect
down-weighting previous observations and can therefore be considered as a
Bayesian analogue of the _forgetting factor_ procedure that has been widely
used in frequentist online literature (see, for example, Haykin, 2002;
Bodenham and Adams, 2017). In this way, the parameter estimates are quicker to
respond to changes in the latent rate structure.

Due to the independence of counts from a Poisson process within non-
overlapping time windows, the posterior density at step $r$ under the BHPP
model factorises as

$\displaystyle p(z,\theta\mid\mathcal{D}_{1:r})$
$\displaystyle\propto\prod_{\ell=1}^{r}p(\mathcal{D}_{\ell}\mid
z,\lambda)\times p(z\mid\pi)\times p(\pi)\times p(\lambda)\propto
p(\mathcal{D}_{r}\mid z,\lambda)\times p(z,\theta\mid\mathcal{D}_{1:(r-1)}),$
(3.5)

where the prior distributions are chosen to be conjugate, as per (2.3) and
(2.5). From (3.5), the posterior distribution at step $r-1$ can be interpreted
as the prior for the posterior distribution at step $r$.

The posterior density in (3.5) is not available in closed form and, due to the
marginal density $p(\mathcal{D}_{1:r})$ being unavailable, can only be
evaluated up to a normalising constant. We thus consider the VB approach with
the mean-field variational family, deploying the CAVI algorithm to compute
approximate posterior distributions. Denoting the approximate posterior at
step $r-1$ by $\hat{q}^{(r-1)}(\theta,z)$, we pass through as the prior for
step $r$ the tempered density

$\displaystyle\hat{q}^{(r-1)}_{\delta}(\theta,z)$
$\displaystyle=\prod_{k,m\in\mathcal{K}}\hat{q}^{(r-1)}_{\delta_{\lambda}}(\lambda_{km})\times\prod_{i\in\mathcal{V}}\hat{q}^{(r-1)}_{\delta_{z}}(z_{i})\times\hat{q}^{(r-1)}_{\delta_{\pi}}(\pi)$
(3.6)
$\displaystyle=\frac{1}{C_{r-1}}\times\prod_{k,m\in\mathcal{K}}\hat{q}^{(r-1)}(\lambda_{km})^{\delta_{\lambda}}\times\prod_{i\in\mathcal{V}}\hat{q}^{(r-1)}(z_{i})^{\delta_{z}}\times\hat{q}^{(r-1)}(\pi)^{\delta_{\pi}},$
(3.7)

controlled by the Bayesian forgetting factor (BFF)
$\delta=(\delta_{\lambda},\delta_{z},\delta_{\pi})\in(0,1]^{3}$. Here,
$C_{r-1}$ is a normalising constant to ensure a valid density. This procedure
is visualised in Figure 4.

Figure 4: A representation of the transformation of the CAVI approximation
from the $r$-th iteration via temperature parameters $\delta$ and its use as a
prior for the $(r+1)$-th iteration. The overlaying of $\mathcal{D}_{r}$ on the
diagonal arrows indicates that the transformed CAVI posterior is combined with
incoming data to yield the next approximate posterior.

The effect of the BFF in down-weighting previous data and placing greater
emphasis on the latest batch $\mathcal{D}_{r}$ is illustrated in Figure 5. A
fully-connected network is simulated on the time interval $[0,5]$ with $N=500$
nodes and $K=2$ groups. An instantaneous rate change in $\lambda_{11}$ is made
at $t=1$. The left panel of Figure 5 shows the results for naive approach
(without a forgetting factor), whereas the right panel includes a BFF and
demonstrates a much faster response to the change and quicker convergence of
the posterior mean to the true value.

Figure 5: Estimated posterior mean of $\lambda\in\mathbb{R}_{+}^{K\times K}$
for a fully connected graph with $N=500$ nodes simulated on $[0,5]$, with
$K=2$ groups, $\pi=(0.6,0.4)$, $\lambda_{22}=8,\ \lambda_{12}=1,\
\lambda_{21}=0.3$, and $\lambda_{11}=2$. At $t=1$, $\lambda_{11}$ changes to
$5$. The four coloured lines are the posterior means of the components
$\lambda$, plotted against iteration. The black dashed horizontal lines are
the true values of $\lambda$ at each time point. $\delta=0.1$ in the right-
hand panel.

###  3.3 Online CAVI updates

We assume that each component of (3.6) takes the same form of its
corresponding complete conditional distribution under the standard BHPP model
in (2.2)–(2.5), raised to a power and normalised:

$\displaystyle\hat{q}^{(r-1)}_{\delta_{\lambda}}(\lambda_{km})$
$\displaystyle=\frac{1}{C_{\lambda_{km},r-1}}\mathrm{Gamma}\left(\alpha_{km}^{(r-1)},\beta_{km}^{(r-1)}\right)^{\delta_{\lambda}},$
$\displaystyle k,m\in\mathcal{K},$
$\displaystyle\hat{q}^{(r-1)}_{\delta_{z}}(z_{i})$
$\displaystyle=\frac{1}{C_{z_{i},r-1}}\mathrm{Categorical}\left(\tau^{(r-1)}\right)^{\delta_{z}},$
$\displaystyle i\in\mathcal{V},$
$\displaystyle\hat{q}^{(r-1)}_{\delta_{\pi}}(\pi)$
$\displaystyle=\frac{1}{C_{\pi,r-1}}\mathrm{Dirichlet}\left(\gamma^{(r-1)}\right)^{\delta_{\pi}}.$

Here, $C_{\lambda_{km},r-1},C_{z_{i},r-1}$ and $C_{\pi,r-1}$ are normalising
constants, ensuring that the densities are valid, and
$\delta=(\delta_{\lambda},\delta_{z},\delta_{\pi})\in(0,1]^{3}$ are
temperature parameters specific to each class of latent variables. Suppose
that data is observed on the interval $I_{r}:=(L_{r-1},L_{r}]$, for
$r\in\mathbb{N}$, and denote the count on edge $(i,j)$ during this interval by
$x_{ij}(I_{r})$. Under the prior structure in (3.6), it is shown in Appendix A
that the CAVI sequential updates are:

  1. 1.

$\hat{q}^{(r)}(\lambda_{km})=\mathrm{Gamma}\left(\alpha^{(r)}_{km},\beta^{(r)}_{km}\right)$
for all $k,m\in\mathcal{K}$ where $\alpha^{(r)}_{km},\beta^{(r)}_{km}$ are
defined as:

$\displaystyle\alpha^{(r)}_{km}$
$\displaystyle=\delta_{\lambda}\left(\alpha^{(r-1)}_{km}-1\right)+\sum_{(i,j)\in\mathcal{E}}\tau_{ik}^{(r-1)}\tau_{jm}^{(r-1)}x_{ij}(I_{r})+1,$
(3.8) $\displaystyle\beta^{(r)}_{km}$
$\displaystyle=\delta_{\lambda}\beta^{(r-1)}_{km}+\Delta\sum_{(i,j)\in\mathcal{E}}\tau_{ik}^{(r-1)}\tau_{jm}^{(r-1)}.$
(3.9)

  2. 2.

$\hat{q}^{(r)}(z_{i})=\mathrm{Categorical}\left(\tau_{i}^{(r)}\right)$ for all
$i\in\mathcal{V}$ where
$\tau_{i}^{(r)}=(\tau_{i1}^{(r)},\dots,\tau_{iK}^{(r)})$, with
$\sum_{k\in\mathcal{K}}\tau_{ik}^{(r)}=1$, satisfies the relation:

$\displaystyle\tau_{ik}^{(r)}$
$\displaystyle\propto\exp\Bigg{\\{}\delta_{z}\left[\psi\left(\gamma^{(r-1)}_{k}\right)-\psi\left(\sum_{m\in\mathcal{K}}\gamma^{(r-1)}_{m}\right)\right]\
+x_{ii}(I_{r})\left[\psi\left(\alpha_{kk}^{(r)}\right)-\log\left(\beta_{kk}^{(r)}\right)\right]-\Delta\frac{\alpha_{kk}^{(r)}}{\beta_{kk}^{(r)}}+$
$\displaystyle\sum_{m\in\mathcal{K}}\Bigg{[}\sum_{j:(i,j)\in\mathcal{E}}\tau_{jm}^{(r)}\bigg{(}x_{ij}(I_{r})\left\\{\psi\left(\alpha^{(r)}_{km}\right)-\log\left(\beta^{(r)}_{km}\right)\right\\}-\Delta\frac{\alpha^{(r)}_{km}}{\beta^{(r)}_{km}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}\
+$
$\displaystyle\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\tau_{j^{\prime}m}^{(r)}\bigg{(}x_{j^{\prime}i}(I_{r})\left\\{\psi\left(\alpha^{(r)}_{mk}\right)-\log\left(\beta^{(r)}_{mk}\right)\right\\}-\Delta\frac{\alpha^{(r)}_{mk}}{\beta^{(r)}_{mk}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j^{\prime})\bigg{)}\Bigg{]}\Bigg{\\}}.$
(3.10)

  3. 3.

$\hat{q}^{(r)}(\pi)=\mathrm{Dirichlet}\left\\{\left(\gamma^{(r)}_{1},\dots,\gamma^{(r)}_{K}\right)\right\\}$
where $\gamma^{(r)}_{k}$ is defined as:

$\displaystyle\gamma^{(r)}_{k}=\delta_{\pi}\left(\gamma_{k}^{(r-1)}-1\right)+\delta_{z}\sum_{i\in\mathcal{V}}\tau_{ik}^{(r)}+1.$
(3.11)

Here $\psi(\cdot)$ is the digamma function. Note that the parameters
$\tau_{ik}^{(r)},\ i\in\mathcal{V},\ k\in\mathcal{K}$, are jointly optimised
via a fixed point solver for (3.10), iteratively moving through the rows of
the matrix $\tau=\\{\tau_{ik}\\}_{i\in\mathcal{V},\
k\in\mathcal{K}}\in[0,1]^{N\times K}$ until convergence. The fixed point
solver is initialised at the $r$-th step by using the value of $\tau$
outputted from the $(r-1)$-th step as the starting point, which also partially
circumvents the problem of label switching (Jasra et al., 2005).

Algorithm 1 describes the full inference procedure based upon the updates in
equations (3.9)–(3.11).

Algorithm 1 Online VB Procedure for BHPP

Initialise $\alpha^{(0)}$, $\beta^{(0)}$, $\gamma^{(0)}$ and set
$\tau_{ik}=1/K$ for all $i\in\mathcal{V}$ and $k\in\mathcal{K}$.

for $r=1,2,\dots$ do

for $\ell=1,2,\dots,N_{\mathrm{CAVI}}$ do

while $\tau$ not converged do

for $j=1$ to $N$ do

Update $\tau_{j}=(\tau_{j1},\dots,\tau_{jK})$ as in equation (3.10),
intialising as the output of update $r-1$.

end for

end while

Normalise $\tau_{i}$ such that $\sum_{k=1}^{K}\tau_{ik}=1$.

Update $\gamma^{(r)}$, $\alpha^{(r)}$ and $\beta^{(r)}$ as in (3.9)–(3.11).

end for

end for

###  3.4 Detecting changepoints

In this work, we aim to detect two types of changepoints: changes in the point
process rates and changes in the group memberships.

To detect changes in the point process rates, we compute the KL-divergence
between the approximate posteriors $\hat{q}^{(r)}(\lambda_{km})$ and
$\hat{q}^{(r-s)}(\lambda_{km})$, for $s=1,2,\dots,\kappa$, where
$\kappa\in\mathbb{N}$ is a lag parameter. We compare the current approximation
with estimates beyond that of the previous update, up to a maximal lag of
$\kappa$, to avoid incorrectly flagging outliers as changes. Figure 6(a)
illustrates the behaviour of the KL-divergence for a lag of 1 and a lag of 2.
In Figure 6(b), we see that using only a lag of 1 would result in the false
identification of a change. There is a trade-off to be made between confidence
and speed of detection. We found that a maximal lag of $\kappa=2$ provided a
good balance.

To compute the KL-divergence, we note that $\hat{q}^{(r)}(\lambda_{km})$ is
gamma distributed for all $r\in\mathbb{N}$. For two random variables
$X_{1}\sim\mathrm{Gamma}(\alpha_{1},\beta_{1})$ and
$X_{2}\sim\mathrm{Gamma}(\alpha_{2},\beta_{2})$ with probability distributions
$p_{1}$ and $p_{2}$ respectively, their KL-divergence is:

$\mathrm{KL}(p_{1}\ \|\
p_{2})=\alpha_{2}\log\frac{\beta_{1}}{\beta_{2}}-\log\frac{\Gamma(\alpha_{1})}{\Gamma(\alpha_{2})}+(\alpha_{1}-\alpha_{2})\psi(\alpha_{1})-(\beta_{1}-\beta_{2})\frac{\alpha_{1}}{\beta_{1}}.$
(3.12)

Two burn-in times are required to make use of the KL-divergence. In
particular, $B_{1}\in\mathbb{N}$ iterations are needed for the algorithm to
converge from initialisation, and $B_{2}\in\mathbb{N}$ further steps needed to
obtain enough samples to flag changes. During the initial $B_{1}+B_{2}$ steps,
stationarity is assumed. After $B_{1}+B_{2}$ steps, we propose to use the
median absolute deviation (MAD) to evaluate the presence of a changepoint, as
it is robust to outliers (Leys et al., 2013). In general, for a dataset
$\mathcal{Y}:=\\{y_{i}\\}_{i=1}^{n}$, the MAD is defined as

$\mathrm{MAD}=\operatorname*{med}\big{(}\\{|y_{i}-\operatorname*{med}(\mathcal{Y})|\\}_{i=1}^{n}\big{)},$
(3.13)

where $\operatorname*{med}(\cdot)$ denotes the median of the elements in a
set. An observation $y_{j}\in\mathcal{Y}$ is classified as an outlier if
$|y_{j}-\operatorname*{med}(\mathcal{Y})|>L\cdot\mathrm{MAD}$, for a choice of
threshold $L\in\mathbb{R}_{+}$.

In our proposed methodology, the elements in the set $\mathcal{Y}$ are the
realised values of the KL-divergence between different gamma distributions,
approximating the posterior distributions for the parameters $\lambda_{km},\
k,m\in\mathcal{K}$ at different time windows. In particular, we write
$y^{r,s}_{km}:=\mathrm{KL}[q^{(r)}(\lambda_{km})\ \|\
q^{(r-s)}(\lambda_{km})]$ and construct
$\mathcal{Y}_{km}^{s}:=\\{y^{r,s}_{km}\\}_{r\geq B_{1}+B_{2}+s}$. A
changepoint is flagged at $t=r\Delta$ if for all $s=1,2,\dots,\kappa$, we
have:

$|y_{km}^{r+s-1,s}-\operatorname*{med}(\mathcal{Y}_{km}^{s})|>L\cdot\mathrm{MAD}.$
(3.14)

It must be remarked that if we choose a maximal lag of $\kappa$, we can only
flag a changepoint a minimum of $\kappa-1$ time steps after the changepoint
has occurred.

We now consider whether the stream must be reset post change. Suppose a change
to the latent rate $\lambda$ occurs at some time $L_{r}<t^{\prime}<L_{r+1}$,
and that pre-change $\lambda\sim\mathrm{Gamma}(\alpha,\beta)$, whereas
$\lambda\sim\mathrm{Gamma}(\alpha^{\prime},\beta^{\prime})$ post-change, for
$\alpha\neq\alpha^{\prime}$ and $\beta\neq\beta^{\prime}$. Pre-change, our
algorithm provides two consecutive CAVI estimates $q^{(r-1)}(\lambda)$ and
$q^{(r)}(\lambda)$, with shapes and rates $\alpha^{(r-1)},\beta^{(r-1)}$ and
$\alpha^{(r)},\beta^{(r)}$, respectively, and post-change, we have CAVI
estimates $q^{(r+1)}(\lambda)$ and $q^{(r+2)}(\lambda)$ with shapes and rates
$\alpha^{(r+1)},\beta^{(r+1)}$ and $\alpha^{(r+2)},\beta^{(r+2)}$. For ease of
analysis, we suppose that
$\alpha^{(r)}-\alpha^{(r-1)}=\alpha^{(r+2)}-\alpha^{(r+1)}=\Delta_{\alpha}$,
and $\beta^{(r)}-\beta^{(r-1)}=\beta^{(r+2)}-\beta^{(r+1)}=\Delta_{\beta}$.

We consider the ratio of the KL-divergences between consecutive estimates pre
and post-change, and expand as $\Delta_{\alpha}\to 0$ to get

$\frac{\mathrm{KL}(q^{(r+2)}\ \|\ q^{(r+1)})}{\mathrm{KL}(q^{(r)}\ \|\
q^{(r-1)})}=\frac{\alpha^{(r+1)}\left(\Delta_{\beta}/\beta^{(r+1)}-\log\left(1+\Delta_{\beta}/\beta^{(r+1)}\right)\right)}{\alpha^{(r-1)}\left(\Delta_{\beta}/\beta^{(r-1)}-\log\left(1+\Delta_{\beta}/\beta^{(r-1)}\right)\right)}+\mathcal{O}(\Delta_{\alpha}).$
(3.15)

Similarly, taking $\Delta_{\beta}\to 0$ yields

$\frac{\mathrm{KL}(q^{(r+2)}\ \|\ q^{(r+1)})}{\mathrm{KL}(q^{(r)}\ \|\
q^{(r-1)})}=\frac{\log\left\\{\Gamma(\alpha^{(r+1)})/\Gamma(\alpha^{(r+1)}+\Delta_{\alpha})\right\\}+\Delta_{\alpha}\psi(\alpha^{(r+1)})}{\log\left\\{\Gamma(\alpha^{(r-1)})/\Gamma(\alpha^{(r-1)}+\Delta_{\alpha})\right\\}+\Delta_{\alpha}\psi(\alpha^{(r-1)})}+\mathcal{O}(\Delta_{\beta}).$
(3.16)

Observation of (3.15) indicates that to first order the ratio depends upon
only the relative size of the difference between their rates, that is
$\Delta_{\beta}/\beta$, and not directly upon the magnitude of $\beta$. This
is not problematic only if the magnitude of $\Delta_{\beta}$ changes linearly
with $\beta$. Also, (3.16) shows that the ratio depends directly on the
magnitude of the scale $\alpha$. It follows that after the identification of a
changepoint at $t=r_{c}\Delta$, $r_{c}\in\mathbb{N}$, the stream of realised
KL-divergence values must be discarded, and the algorithm wait a further
$B_{2}$ steps to obtain enough samples for the new set
$\mathcal{Y}_{km}^{s}:=\\{y_{km}^{r,s}\\}_{r\geq r_{c}+B_{2}+s}$ which is used
to detect any subsequent changes. A further benefit of consider a maximal lag
$\kappa>1$ is that more KL-divergence samples are obtained for the same number
of observations.

To detect changes in the group memberships, the approximation
$\prod_{i\in\mathcal{V}}\hat{q}(z_{i})$ to the posterior of the latent group
memberships in (3.10) provides a node-level vector $\tau_{i}$ ($i=1,...,N)$ of
group membership probabilities. Computing
$\operatorname*{argmax}_{k}\tau_{ik}$ provides a group assignment for node
$i$, and by comparing assignments between runs, one has a natural way of
flagging changes. However, this approach is unsatisfactory as it does not
account for the magnitude of the change in probability. The KL-divergence also
cannot be used for flagging group changes, since $\mathrm{KL}(p_{1}\ \|\
p_{2})$ is only defined only when $p_{2}>0$ for all values in the support of
$p_{1}$. Instead, we utilise the Jensen-Shannon divergence (JS; Lin, 1991),
defined for two distributions $p_{1},p_{2}$ as

$\mathrm{JS}(p_{1}\ \|\ p_{2})=\frac{1}{2}\mathrm{KL}(p_{1}\ \|\
p)+\frac{1}{2}\mathrm{KL}(p_{2}\ \|\ p),$ (3.17)

where $p:=(p_{1}+p_{2})/2$ is a 50-50 mixture distribution of $p_{1}$ and
$p_{2}$. The JS-divergence is symmetric in its arguments, and avoids the
undefined values encountered with the KL-divergence. For
$q_{1},q_{2}\in\Delta^{n}$, where $\Delta^{n}$ denotes the $n$-dimensional
simplex, the JS-divergence between discrete random variables
$X_{1}\sim\mathrm{Categorical}(q_{1})$ and
$X_{2}\sim\mathrm{Categorical}(q_{2})$, with probability mass functions
$p_{1}$ and $p_{2}$ respectively, takes the following form:

$\mathrm{JS}(p_{1}\ \|\
p_{2})=\frac{1}{2}\left(\sum_{i=1}^{n}q_{1i}\log\frac{q_{1i}}{(q_{1i}+q_{2i})/2}+\sum_{i=1}^{n}q_{2i}\log\frac{q_{2i}}{(q_{1i}+q_{2i})/2}\right).$
(3.18)

The MAD is then used to detect changes to the stream of logged values, but
with two additional requirements on the group probabilities:

$\displaystyle\operatorname*{argmax}_{k}\tau^{(r)}_{ik}\neq\operatorname*{argmax}_{k}\tau^{(r-s)}_{ik},$
$\displaystyle\text{for all }s=1,\dots,\kappa,$ (3.19)
$\displaystyle\operatorname*{argmax}_{k}\tau^{(r-s_{1})}_{ik}=\operatorname*{argmax}_{k}\tau^{(r-s_{2})}_{ik},$
$\displaystyle\text{for all }s_{1},s_{2}\in\\{1,\dots,\kappa\\}.$ (3.20)

If the MAD condition and both (3.19) and (3.20) are met, then a changepoint is
flagged. The conditions (3.19) and (3.20) ensure that changes are flagged only
if the most likely group assignment changes after a window of $\kappa$ steps
where the most likely group assignment was stable. It should noted that for
this stream, as the probabilities are constrained to sum to 1, we do not reset
these values after a changepoint is detected.

(a) A true changepoint. (b) An outlier. Figure 6: Illustration of the
behaviour of the KL-divergence with different lags for a changepoint and an
outlier.

###  3.5 Online VB for the dynamic BHPP with an unknown underlying graph

Until now, it has been assumed that the graph is fully connected. We relax
this assumption and extend the framework to the setting where
$\bm{A}=\\{a_{ij}\\}\in\\{0,1\\}^{N\times N}$ is unknown. Assuming $a_{ij}=1$
when there is no edge adversely affects estimation of the latent rates;
instead of the absence of an event indicating the absence of an edge, it
instead just leads to a lower rate estimate. This motivates a method that can
take into account the possibility that $a_{ij}=0$. The work here builds upon
the sparse setup of Matias et al. (2018).

For each $(i,j)\in\mathcal{R}$, an additional latent variable
$a_{ij}\in\\{0,1\\}$ is introduced. A new latent group membership
$z_{i}^{\prime}\in\mathcal{K}^{\prime}:=\\{1,\dots,K^{\prime}\\}$ is assigned
to each $i\in\mathcal{V}$, and $a_{ij}$ is Bernoulli-distributed conditional
upon $z_{i}^{\prime}$ and $z_{j}^{\prime}$, according to a stochastic
blockmodel (SBM; Holland et al., 1983). Each $a_{ij}$ captures the probability
that $(i,j)\in\mathcal{R}$ is also an element of $\mathcal{E}$, and functions
to reduce the contribution of edges with no arrivals to the estimates of the
underlying rates. The new model is then expressed as:

$\displaystyle x_{ij}(t)\mid a_{ij},z_{i},z_{j},\lambda_{z_{i}z_{j}}$
$\displaystyle\sim\mathrm{Poisson}(a_{ij}\lambda_{z_{i}z_{j}}t),$
$\displaystyle\text{for all }i,j\in\mathcal{V},$ (3.21)
$\displaystyle\lambda_{km}$
$\displaystyle\sim\mathrm{Gamma}\left(\alpha_{km}^{0},\beta_{km}^{0}\right),$
$\displaystyle\text{for all }k,m\in\mathcal{K},$ (3.22) $\displaystyle
a_{ij}\mid z_{i}^{\prime},z_{j}^{\prime},\rho_{z_{i}^{\prime}z_{j}^{\prime}}$
$\displaystyle\sim\mathrm{Bernoulli}(\rho_{z_{i}^{\prime}z_{j}^{\prime}}),$
$\displaystyle\text{for all }i,j\in\mathcal{V},$ (3.23)
$\displaystyle\rho_{k^{\prime}m^{\prime}}$
$\displaystyle\sim\mathrm{Beta}\left(\eta_{k^{\prime}m^{\prime}}^{0},\zeta_{k^{\prime}m^{\prime}}^{0}\right),$
$\displaystyle\text{for all }k^{\prime},m^{\prime}\in\mathcal{K}^{\prime},$
(3.24) $\displaystyle z_{i}\mid\pi$
$\displaystyle\sim\mathrm{Categorical}(\pi),$ $\displaystyle\text{for all
}i\in\mathcal{V},$ (3.25) $\displaystyle z_{i}^{\prime}\mid\mu$
$\displaystyle\sim\mathrm{Categorical}(\mu),$ $\displaystyle\text{for all
}i\in\mathcal{V},$ (3.26) $\displaystyle\pi$
$\displaystyle\sim\mathrm{Dirichlet}\left(\gamma^{0}\right),$ (3.27)
$\displaystyle\mu$ $\displaystyle\sim\mathrm{Dirichlet}\left(\xi^{0}\right),$
(3.28)

We allow the most general setting in which the group structure that governs
the rates is separate from that which determines edge connection
probabilities. We assume that the graph is static, and so changes to the
memberships that drive the edge-processes should not affect the determination
of the graph adjacency matrix. The decoupling of the edge-process and edge-
connection group memberships ensure this. Furthermore, as has been noted, SBMs
provide good approximations to any exchangeable random graph model (Airoldi et
al., 2013), and they therefore offer a logical model for the graph generation
mechanism. This new model structure is represented in graphical model form in
Figure 7.

Figure 7: A directed acyclic graphical model representation of the model in
(3.21)-(3.28).

Given the graph is considered to be static, if $x_{ij}(t)>0$ for some $t$,
then $(i,j)\in\mathcal{E}$ (corresponding to $a_{ij}=1$) with probability one.
Otherwise, if $x_{ij}(t)=0$ then we either have no edge (corresponding to
$a_{i}j=0$), or a non-null intensity process with no events in $(0,t]$. The
difference with the Section 2.1 is that we now have the additional local
latent variables, $z$, $z^{\prime}$ and $\bm{A}$, all of which grow with $N$
(linearly for $z$ and $z^{\prime}$, and quadratically for $\bm{A}$).

To implement variational inference, the full conditional distribution
$p(a_{ij}\mid
z_{i},z_{j},z_{i}^{\prime},z_{j}^{\prime},\lambda_{z_{i}z_{j}},\rho_{z_{i}^{\prime}z_{j}^{\prime}},\mathcal{D}_{1:r})$
must be written down exactly, as per (3.4). Using the fact that inter-arrival
times from a Poisson process are exponentially distributed, Matias et al.
(2018) shows that the conditional posterior probability of an edge
$(i,j)\in\mathcal{E}$ in the case of static $\lambda$ takes the form

$p(a_{ij}=1\mid\mathcal{D}_{1:r},z_{i}=k,z_{j}=m,z_{i}^{\prime}=k^{\prime},z_{j}^{\prime}=m^{\prime},\lambda_{km},\rho_{k^{\prime}m^{\prime}})=\\\
\mathbb{I}\\{x_{ij}(r\Delta)>0\\}+\frac{\mathbb{I}\\{x_{ij}(r\Delta)=0\\}\rho_{k^{\prime}m^{\prime}}\exp\left(-\lambda_{km}r\Delta\right)}{1-\rho_{k^{\prime}m^{\prime}}+\rho_{k^{\prime}m^{\prime}}\exp\left(-\lambda_{km}r\Delta\right)}.$
(3.29)

Since $\bm{A}$ is static, the posterior for each $a_{ij}$ should be computed
using the full data $\mathcal{D}_{1:r}$. Note that the conditional
distribution in (3.29) depends on $\mathcal{D}_{1:r}$ only via the counting
process $x_{ij}(\cdot)$ at time $r\Delta$.

For the model in (3.21)-(3.28), the true posterior distribution after $r$
batches factorises as

$\displaystyle p(\bm{A},z,z^{\prime},\theta\mid\mathcal{D}_{1:r})$
$\displaystyle\propto p(\mathcal{D}_{r}\mid\bm{A},z,z^{\prime}\theta)\times
p(\bm{A},z,z^{\prime},\theta\mid\mathcal{D}_{1:(r-1)}),$ (3.30)

similarly to the posterior distribution in (3.5). Therefore, inference on this
model can be performed in much the same way as inference on the original
model. However, as the true conditional posterior for $\bm{A}$ is tractable,
the update procedure is slightly different. Unlike Matias et al. (2018), we
cannot simply take $q^{(r)}(a_{ij})$ to be Bernoulli distributed with
probability as in (3.29) as $\lambda$ is, in general, not static. Instead, to
account for variations in $\lambda$ over the full observation window, we
propose an approximation to the posterior distribution of $a_{ij}$ at
$t=r\Delta$, conditional upon
$z_{i}=k,z_{j}=m,z_{i}^{\prime}=k^{\prime},z_{j}^{\prime}=m^{\prime}$, of the
form $\hat{q}^{(r)}(a_{ij})=\mathrm{Bernoulli}\left(\sigma_{ij}^{(r)}\right)$,
where

$\sigma_{ij}^{(r)}=\mathbb{I}\\{x_{ij}(r\Delta)>0\\}+\frac{\mathbb{I}\\{x_{ij}(r\Delta)=0\\}\hat{\rho}_{k^{\prime}m^{\prime}}^{(r)}\exp\left(-\Delta\sum_{\ell=0}^{r}\hat{\lambda}^{(\ell)}_{km}\right)}{1-\hat{\rho}_{k^{\prime}m^{\prime}}^{(r)}+\hat{\rho}_{k^{\prime}m^{\prime}}^{(r)}\exp\left(-\Delta\sum_{\ell=0}^{r}\hat{\lambda}^{(\ell)}_{km}\right)}.$
(3.31)

Here, $\hat{\rho}_{k^{\prime}m^{\prime}}^{(\ell)}$ and
$\hat{\lambda}_{km}^{(\ell)}$ denote the mean of our approximate posterior
distribution for $\rho_{k^{\prime}m^{\prime}}^{(\ell)}$ and
$\lambda_{km}^{(\ell)}$, respectively, at update $\ell$. At update $r$, we
propose a variational approximation of the form
$q^{(r)}(\bm{A},z,z^{\prime},\theta)=\prod_{(i,j)\in\mathcal{R}}\mathrm{Bernoulli}(\sigma_{ij}^{(r)})\times
q(\theta,z,z^{\prime})$, where $\sigma_{ij}^{(r)}$ is as in (3.31) and
$q(\theta,z,z^{\prime})\in\mathcal{F}^{2N+S}$. A two-step estimation procedure
follows naturally, wherein $q(\theta,z,z^{\prime})$ is approximated using CAVI
as was done previously, and then $\sigma^{(r)}_{ij}$ is updated using (3.31)
and the CAVI approximations to the posteriors of $\lambda$ and $\rho$. Note
that since $a_{ij}$ for $i,j\in\mathcal{V}$, $z_{i}^{\prime}$ for
$i\in\mathcal{V}$, $\rho_{k^{\prime}m^{\prime}}$ for
$k^{\prime},m^{\prime}\in\mathcal{K}^{\prime}$ and $\mu$ are considered static
parameters, forgetting factors for these quantities are not required.

Consider the BHPP model with unknown graph structure, as in (3.21)-(3.28),
with global parameters $\theta=(\lambda,\pi,\mu,\rho)$ and local parameters
$z$, $z^{\prime}$ and $\bm{A}$. At step $r-1$, we approximate the posterior
density $p(\bm{A},z,z^{\prime},\theta\mid\mathcal{D}_{1:(r-1)})$ by
$\hat{q}^{(r-1)}(\bm{A},z,z^{\prime},\theta)$, which is the product of the
optimal CAVI solution $\hat{q}^{(r-1)}(z,z^{\prime},\theta)$ and the fixed-
form approximation $\hat{q}^{(r-1)}(\bm{A})$ of (3.31). For the BHPP model
with known graph structure, all latent variables could change, and thus in
(3.7) all components were tempered before being passed as the prior for step
$r$. However, in the case of an unknown graph structure, as the graph is
assumed to be static, only the components of the dynamic latent variables are
tempered when constructing the prior for step $r$. We pass through as the
prior for step $r$ the partially tempered density

$\displaystyle\hat{q}^{(r-1)}_{\delta}(\bm{A},z,\theta)=$
$\displaystyle\prod_{k,m\in\mathcal{K}}\hat{q}^{(r-1)}(\rho_{km})\times\prod_{i\in\mathcal{V}}\hat{q}^{(r-1)}(z_{i}^{\prime})\times\prod_{(i,j)\in\mathcal{R}}\hat{q}^{(r-1)}(a_{ij})\times\hat{q}^{(r-1)}(\mu)\times$
(3.32)
$\displaystyle\prod_{i\in\mathcal{V}}\hat{q}^{(r-1)}_{\delta_{z}}(z_{i})\times\prod_{k,m\in\mathcal{K}}\hat{q}^{(r-1)}_{\delta_{\lambda}}(\lambda_{km})\times\hat{q}^{(r-1)}_{\delta_{\pi}}(\pi).$
(3.33)

We assume that each component of (3.33) takes the same form as its
corresponding complete conditional distribution under the BHPP model with
unknown graph structure of (3.21)-(3.28). Additionally, the distributions for
$\lambda$, $\pi$ and $z$ are raised to a power and normalised. In summary:

$\displaystyle\hat{q}^{(r-1)}(a_{ij})$
$\displaystyle=\mathrm{Bernoulli}\left(\sigma_{ij}^{{(r-1)}}\right),$
$\displaystyle i,j\in\mathcal{V},$
$\displaystyle\hat{q}^{(r-1)}_{\delta_{\lambda}}(\lambda_{km})$
$\displaystyle=\frac{1}{C_{\lambda_{km},r-1}}\mathrm{Gamma}\left(\alpha_{km}^{(r-1)},\beta_{km}^{(r-1)}\right)^{\delta_{\lambda}},$
$\displaystyle k,m\in\mathcal{K},$
$\displaystyle\hat{q}^{(r-1)}(\rho_{k^{\prime}m^{\prime}})$
$\displaystyle=\mathrm{Beta}\left(\eta^{(r-1)}_{k^{\prime}m^{\prime}},\zeta^{(r-1)}_{k^{\prime}m^{\prime}}\right),$
$\displaystyle k^{\prime},m^{\prime}\in\mathcal{K}^{\prime},$
$\displaystyle\hat{q}^{(r-1)}(z_{i})$
$\displaystyle=\frac{1}{C_{z_{i},r-1}}\mathrm{Categorical}\left(\tau_{i}^{(r-1)}\right)^{\delta_{z}},$
$\displaystyle i\in\mathcal{V},$
$\displaystyle\hat{q}^{(r-1)}(z_{i}^{\prime})$
$\displaystyle=\mathrm{Categorical}\left(\nu_{i}^{(r-1)}\right),$
$\displaystyle i\in\mathcal{V},$
$\displaystyle\hat{q}^{(r-1)}_{\delta_{\pi}}(\pi)$
$\displaystyle=\frac{1}{C_{\pi,r-1}}\mathrm{Dirichlet}\left(\gamma^{(r-1)}\right)^{\delta_{\pi}},$
$\displaystyle\hat{q}^{(r-1)}(\mu)$
$\displaystyle=\mathrm{Dirichlet}\left(\xi^{(r-1)}\right).$

Suppose that data is again observed on the interval $I_{r}:=(L_{r-1},L_{r}]$,
for $r\in\mathbb{N}$. Under the prior structure in (3.33), the CAVI sequential
updates take the form:

  1. 1.

$\hat{q}^{(r)}(\lambda_{km})=\mathrm{Gamma}\left(\alpha_{km}^{(r)},\beta_{km}^{(r)}\right)$
where for all $k,m\in\mathcal{K}$ we define $\alpha_{km}^{(r)}$ and
$\beta_{km}^{(r)}$ as

$\displaystyle\alpha_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\left(\alpha^{(r-1)}_{km}-1\right)+\sum_{(i,j)\in\mathcal{R}}\tau^{(r-1)}_{ik}\tau^{(r-1)}_{jm}\sigma_{ij}^{(r-1)}x_{ij}^{(r)}+1,$
(3.34) $\displaystyle\beta_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\beta^{(r-1)}_{km}+\sum_{(i,j)\in\mathcal{R}}\tau^{(r-1)}_{ik}\tau^{(r-1)}_{jm}\sigma_{ij}^{(r-1)}.$
(3.35)

  2. 2.

$\hat{q}^{(r)}(\rho_{k^{\prime}m^{\prime}})=\mathrm{Beta}\left(\eta^{(r)}_{k^{\prime}m^{\prime}},\zeta^{(r)}_{k^{\prime}m^{\prime}}\right)$
where for all $k^{\prime},m^{\prime}\in\mathcal{K}^{\prime}$, we define
$\eta_{k^{\prime}m^{\prime}}^{(r)}$ and $\zeta_{k^{\prime}m^{\prime}}^{(r)}$
as

$\displaystyle\eta_{k^{\prime}m^{\prime}}^{(r)}$
$\displaystyle=\eta_{k^{\prime}m^{\prime}}^{(r-1)}+\sum_{(i,j)\in\mathcal{R}}\nu^{(r-1)}_{ik^{\prime}}\nu^{(r-1)}_{jm^{\prime}}\sigma_{ij}^{(r-1)}+1,$
(3.36) $\displaystyle\zeta_{k^{\prime}m^{\prime}}^{(r)}$
$\displaystyle=\zeta_{k^{\prime}m^{\prime}}^{(r-1)}+\sum_{(i,j)\in\mathcal{R}}\nu^{(r-1)}_{ik^{\prime}}\nu^{(r-1)}_{jm^{\prime}}\left(1-\sigma_{ij}^{(r-1)}\right).$

  3. 3.

$\hat{q}^{(r)}(z_{i})=\mathrm{Categorical}\left(\tau_{i}^{(r)}\right)$ for all
$i\in\mathcal{V}$ where
$\tau_{i}^{(r)}=(\tau_{i1}^{(r)},\dots,\tau_{iK}^{(r)})$, with
$\sum_{k\in\mathcal{K}}\tau_{ik}^{(r)}=1$, satisfies the relation:

$\displaystyle\tau_{ik}^{(r)}$
$\displaystyle\propto\exp\Bigg{\\{}\delta_{z}\left[\psi\left(\gamma^{(r-1)}_{k}\right)-\psi\left(\sum_{\ell=1}^{K}\gamma^{(r-1)}_{\ell}\right)\right]+\sum_{j\in\mathcal{V}}\sum_{m\in\mathcal{K}}\tau_{jm}^{(r)}\Bigg{[}-\Delta\left(\sigma_{ij}^{(r-1)}\frac{\alpha^{(r)}_{km}}{\beta^{(r)}_{km}}+\sigma_{ji}^{(r-1)}\frac{\alpha^{(r)}_{mk}}{\beta^{(r)}_{mk}}\right)+$
$\displaystyle\qquad x_{ij}(I_{r})\
\sigma_{ij}^{(r-1)}\left[\psi\left(\alpha^{(r)}_{km}\right)-\log\left(\beta^{(r)}_{km}\right)\right]+x_{ji}(I_{r})\
\sigma_{ji}^{(r-1)}\left[\psi\left(\alpha^{(r)}_{mk}\right)-\log\left(\beta^{(r)}_{mk}\right)\right]\Bigg{]}\times$
(3.37)
$\displaystyle\qquad\bigg{(}1-\mathbb{I}_{\\{i\\}}(j)\mathbb{I}_{\\{k\\}}(m)\bigg{)}+\sigma_{ii}^{(r-1)}\bigg{[}x_{ii}(I_{r})\left[\psi\left(\alpha^{(r)}_{kk}\right)-\log\left(\beta^{(r)}_{kk}\right)\right]-\Delta\frac{\alpha_{kk}}{\beta_{kk}}\bigg{]}\Bigg{\\}}.$
(3.38)

  4. 4.

$\hat{q}^{(r)}(z_{i}^{\prime})=\mathrm{Categorical}\left(\nu_{i}^{(r)}\right)$
for all $i\in\mathcal{V}$ where
$\nu_{i}^{(r)}=(\nu_{i1}^{(r)},\dots,\nu_{iK^{\prime}}^{(r)})$, with
$\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\nu_{ik^{\prime}}^{(r)}=1$, satisfies
the relation:

$\displaystyle\nu_{ik^{\prime}}^{(r)}$
$\displaystyle\propto\exp\Bigg{\\{}\psi\left(\xi^{(r)}_{k^{\prime}}\right)-\psi\left(\sum_{\ell^{\prime}=1}^{K^{\prime}}\xi^{(r)}_{\ell^{\prime}}\right)+\sum_{j\in\mathcal{V}}\sum_{m^{\prime}\in\mathcal{K}^{\prime}}\nu_{jm^{\prime}}^{(r)}\Bigg{[}\sigma_{ij}^{(r)}\left[\psi\left(\eta^{(r)}_{k^{\prime}m^{\prime}}\right)-\psi\left(\eta^{(r)}_{k^{\prime}m^{\prime}}+\zeta^{(r)}_{k^{\prime}m^{\prime}}\right)\right]+$
$\displaystyle\qquad\left(1-\sigma_{ij}^{(r)}\right)\left[\psi\left(\zeta^{(r)}_{k^{\prime}m^{\prime}}\right)-\psi\left(\zeta_{k^{\prime}m^{\prime}}^{(r)}+\eta_{k^{\prime}m^{\prime}}^{(r)}\right)\right]+\sigma_{ji}^{(r)}\left[\psi\left(\eta^{(r)}_{m^{\prime}k^{\prime}}\right)-\psi\left(\eta^{(r)}_{m^{\prime}k^{\prime}}+\zeta^{(r)}_{m^{\prime}k^{\prime}}\right)\right]+$
$\displaystyle\qquad\left(1-\sigma_{ji}^{(r)}\right)\left[\psi\left(\zeta^{(r)}_{m^{\prime}k^{\prime}}\right)-\psi\left(\zeta_{m^{\prime}k^{\prime}}^{(r)}+\eta_{m^{\prime}k^{\prime}}^{(r)}\right)\right]\Bigg{]}\bigg{(}1-\mathbb{I}_{\\{k^{\prime}\\}}(m^{\prime})\mathbb{I}_{\\{i\\}}(j)\bigg{)}+$
(3.39)
$\displaystyle\qquad\sigma_{ii}^{(r)}\bigg{[}\psi\left(\eta^{(r)}_{k^{\prime}k^{\prime}}\right)-\psi\left(\eta^{(r)}_{k^{\prime}k^{\prime}}+\zeta^{(r)}_{k^{\prime}k^{\prime}}\right)\bigg{]}+\left(1-\sigma_{ii}^{(r)}\right)\bigg{[}\psi\left(\zeta^{(r)}_{k^{\prime}k^{\prime}}\right)-\psi\left(\eta^{(r)}_{k^{\prime}k^{\prime}}+\zeta^{(r)}_{k^{\prime}k^{\prime}}\right)\bigg{]}\Bigg{\\}}.$
(3.40)

  5. 5.

$\hat{q}^{(r)}(\pi)=\mathrm{Dirichlet}\left(\gamma^{(r)}\right)$ where for all
$k\in\mathcal{K}$, we define $\gamma^{(r)}_{k}$ as:

$\displaystyle\gamma^{(r)}_{k}=\delta_{\pi}\left(\gamma_{k}^{(r-1)}-1\right)+\delta_{z}\sum_{i\in\mathcal{V}}\tau^{(r)}_{ik}+1,\quad
k\in\mathcal{K}.$ (3.41)

  6. 6.

$\hat{q}^{(r)}(\mu)=\mathrm{Dirichlet}\left(\xi^{(r)}\right)$ where for all
$k^{\prime}\in\mathcal{K}^{\prime}$, we define $\xi_{k^{\prime}}^{(r)}$ as:

$\displaystyle\xi^{(r)}_{k^{\prime}}=\xi_{k^{\prime}}^{(r-1)}+\sum_{i\in\mathcal{V}}\nu^{(r)}_{ik^{\prime}},\quad
k^{\prime}\in\mathcal{K}^{\prime}.$ (3.42)

  7. 7.

$\hat{q}^{(r)}(a_{ij})=\mathrm{Bernoulli}\left(\sigma_{ij}^{(r)}\right)$,
where $\sigma_{ij}^{(r)}$ is defined in (3.31).

###  3.6 Online VB for the dynamic BHPP with an unknown number of groups

In real-world applications, the number of groups $K$ is usually unknown, and
it must be estimated from the observed data. To address this issue, we adopt a
Bayesian non-parametric approach, proposing a Dirichlet process prior
(Ferguson, 1973) on the group memberships. In particular, we replace the
Dirichlet prior distributions in (2.5) and (3.27) with a Griffiths-Engen-
McCloskey (GEM; Pitman, 2002) prior distribution, with parameter $\nu$,
written $\pi\sim\mathrm{GEM}(\nu)$. This prior distribution corresponds to an
infinite limit of a Dirichlet distribution:
$\mathrm{GEM}(\nu)=\lim_{K\to\infty}\mathrm{Dirichlet}(\nu\bm{1}_{K}/K)$,
where $\bm{1}_{K}$ is the vector of ones of length $K$. The full model
becomes:

$\displaystyle x_{ij}(t)\mid z_{i},z_{j},\lambda_{z_{i}z_{j}}$
$\displaystyle\sim\mathrm{Poisson}(\lambda_{z_{i}z_{j}}t),$
$\displaystyle\text{for all }(i,j)\in\mathcal{E},$ (3.43)
$\displaystyle\lambda_{km}$ $\displaystyle\sim\mathrm{Gamma}(\alpha,\beta),$
$\displaystyle\text{for all }k,m=1,2,\dots,$ (3.44) $\displaystyle
z_{i}\mid\pi$ $\displaystyle\sim\mathrm{Categorical}(\pi),$
$\displaystyle\text{for all }i\in\mathcal{V},$ (3.45) $\displaystyle\pi$
$\displaystyle\sim\mathrm{GEM}(\nu),$ (3.46)

where $\pi$ represents an infinite sequence $\pi_{1},\pi_{2},\dots$ such that
$\pi_{k}\geq 0$ for all $k=1,2,\dots$ and $\sum_{k=1}^{\infty}\pi_{k}=1$. The
GEM prior distribution also corresponds to the distribution of proportions
obtained under a stick-breaking representation (Sethuraman, 1994) of a
Dirichlet process. Therefore, the proportions $\pi$ can be reparametrised as a
product of variables $u_{1},u_{2},\dots\in[0,1]$ drawn from independent beta
distributions, as follows:

$\displaystyle u_{k}\sim\mathrm{Beta}(1,\nu),$
$\displaystyle\pi_{k}:=u_{k}\prod_{\ell=1}^{k-1}(1-u_{\ell}),$
$\displaystyle\text{for }k=1,2,\dots.$ (3.47)

This decomposition is particularly useful to derive an online variational
inference algorithm for the BHPP model with GEM priors (GEM-BHPP), following
Blei and Jordan (2006). In particular, within a mean-field approximation
$q(\lambda,u,z)=q(\lambda)\times q(u)\times q(z)$ for the posterior
distribution $p(\lambda,u,z\mid\mathcal{D}_{1:r})$, a variational
approximation $q(u)=\prod_{\ell=1}^{\infty}q(u_{\ell})$ is posited directly on
$u_{1},u_{2},\dots$, rather than $\pi$. This approximation is truncated at a
level $L\in\mathbb{N}$, implying that $q(u_{\ell})=\delta_{0}(u_{\ell})$ for
$\ell=L+1,L+2,\dots$, effectively resulting in an $L$-dimensional probability
vector from (3.47). It should be pointed out how this approach differs from
the model set-up we considered earlier. In the GEM-BHPP, we make no truncation
in the model set-up, but rather only in the variational approximation to it.
This is in contrast to the method we presented in Sections 2.1, wherein one
would be using a finite-dimensional Dirichlet distribution directly in the
prior structure.

We introduce the notation $\omega_{i}^{(r)}$, $\nu_{i}^{(r)}$,
$\alpha_{km}^{(r)}$ and $\beta_{km}^{(r)}$, for $i,k,m\in\\{1,\dots,L\\}$. If
we define $\omega^{(0)}_{i}\equiv 1$, $\nu_{i}^{(0)}\equiv\nu$ for all
$i\in\\{1,\dots,L\\}$ and $\alpha_{km}^{(0)}\equiv\alpha$ and
$\beta_{km}^{(0)}\equiv\beta$, for all $k,m\in\\{1,\dots,L\\}$, the first
iteration takes the form of the model specified in (3.43)–(3.46). Using this
notation, we can formulate our inference procedure based on the same process
of sequential CAVI updates.

Consider the GEM-BHPP model in (3.46)–(3.43) with a stick-breaking
reparametrisation of $\pi$ as in 3.47, global parameters $\theta=(\lambda,u)$
and local parameters $z$. The posterior distribution
$p(\theta,z\mid\mathcal{D}_{1:(r-1)})$ is approximated by the optimal CAVI
solution $\hat{q}^{(r-1)}(\theta,z)$. Similarly to (3.7), this CAVI
approximation is tempered to form $\hat{q}^{(r-1)}_{\delta}(\theta,z)$, taking
the following product form:

$\hat{q}^{(r-1)}_{\delta}(\theta,z)=\prod_{k=1}^{L}\prod_{m=1}^{L}\hat{q}^{(r-1)}_{\delta_{\lambda}}(\lambda_{km})\times\prod_{i\in\mathcal{V}}\hat{q}^{(r-1)}_{\delta_{z}}(z_{i})\times\prod_{\ell=1}^{L}\hat{q}^{(r-1)}_{\delta_{u}}(u_{\ell}).$
(3.48)

This tempered density is then passed through as the prior for step $r$. We
assume that each component of (3.48) takes the same form of its corresponding
complete conditional distribution under the GEM-BHPP model in (3.43)–(3.46),
raised to a power and normalised:

$\displaystyle\hat{q}^{(r-1)}_{\delta_{\lambda}}(\lambda_{km})$
$\displaystyle=\frac{1}{C_{\lambda_{km},r-1}}\mathrm{Gamma}\left(\alpha_{km}^{(r-1)},\beta_{km}^{(r-1)}\right)^{\delta_{\lambda}},$
$\displaystyle k,m\in\\{1,\dots,L\\},$ (3.49)
$\displaystyle\hat{q}^{(r-1)}_{\delta_{z}}(z_{i})$
$\displaystyle=\frac{1}{C_{z_{i},r-1}}\mathrm{Categorical}(\tau^{(r-1)})^{\delta_{z}},$
$\displaystyle i\in\mathcal{V},$ (3.50)
$\displaystyle\hat{q}^{(r-1)}_{\delta_{u}}(u_{\ell})$
$\displaystyle=\frac{1}{C_{u_{\ell},r-1}}\mathrm{Beta}\left(\omega_{\ell}^{(r-1)},\nu_{\ell}^{(r-1)}\right)^{\delta_{u}},$
$\displaystyle\ell\in\\{1,\dots,L\\},$ (3.51)

where $C_{\lambda_{km},r-1},C_{z_{i},r-1}$ and $C_{u_{\ell},r-1}$ are
normalising constants, ensuring that the densities are valid, and
$\delta=(\delta_{\lambda},\delta_{z},\delta_{u})\in(0,1]^{3}$ are temperature
parameters specific to each class of latent variables. Suppose that data is
observed on the interval $I_{r}:=(L_{r-1},L_{r}]$, for $r\in\mathbb{N}$, and
denote the count on edge $(i,j)$ during this interval by $x_{ij}(I_{r})$.
Under the prior structure in (3.48), the CAVI sequential updates then take the
form

  1. 1.

$\hat{q}^{(r)}(\lambda_{km})=\mathrm{Gamma}(\alpha^{(r)}_{km},\beta^{(r)}_{km})$
for all $k,m\in\mathcal{K}$ where $\alpha^{(r)}_{km},\beta^{(r)}_{km}$ are
defined as:

$\displaystyle\alpha_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\left(\alpha_{km}^{(r-1)}-1\right)+\sum_{(i,j)\in\mathcal{E}}\tau_{ik}^{(r-1)}\tau_{jm}^{(r-1)}x_{ij}(I_{r})+1,$
(3.52) $\displaystyle\beta_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\beta_{km}^{(r-1)}+\Delta\sum_{(i,j)\in\mathcal{E}}\tau_{ik}^{(r-1)}\tau_{jm}^{(r-1)}.$
(3.53)

  2. 2.

$\hat{q}^{(r)}(z_{i})=\mathrm{Categorical}\left(\tau_{i}^{(r)}\right)$ for all
$i\in\mathcal{V}$ where
$\tau_{i}^{(r)}=(\tau_{i1}^{(r)},\dots,\tau_{iL}^{(r)})$, with
$\sum_{k\in\mathcal{L}}\tau_{ik}^{(r)}=1$, satisfies the relation:

$\displaystyle\tau_{ik}^{(r)}\propto\exp$
$\displaystyle\Bigg{\\{}\delta_{z}\left[\psi\left(\omega_{k}^{(r-1)}\right)-\psi\left(\omega_{k}^{(r-1)}+\nu_{k}^{(r-1)}\right)+\sum_{\ell=1}^{k-1}\bigg{(}\psi\left(\nu_{\ell}^{(r-1)}\right)-\psi\left(\omega_{\ell}^{(r-1)}+\nu_{\ell}^{(r-1)}\right)\bigg{)}\right]+$
$\displaystyle\sum_{m=1}^{L}\Bigg{[}\sum_{j:(i,j)\in\mathcal{E}}\tau_{jm}^{(r)}\bigg{(}x_{ij}(I_{r})\left[\psi\left(\alpha_{km}^{(r)}\right)-\log\left(\beta_{km}^{(r)}\right)\right]-\Delta\frac{\alpha_{km}^{(r)}}{\beta_{km}^{(r)}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}\
+$
$\displaystyle\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\tau_{j^{\prime}m}^{(r)}\bigg{(}x_{j^{\prime}i}(I_{r})\left[\psi\left(\alpha_{mk}^{(r)}\right)-\log\left(\beta_{mk}^{(r)}\right)\right]-\Delta\frac{\alpha_{mk}^{(r)}}{\beta_{mk}^{(r)}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j^{\prime})\bigg{)}\Bigg{]}+$
(3.54) $\displaystyle
x_{ii}(I_{r})\left[\psi\left(\alpha_{kk}^{(r)}\right)-\log\left(\beta_{kk}^{(r)}\right)\right]-\Delta\frac{\alpha_{kk}^{(r)}}{\beta_{kk}^{(r)}}\Bigg{\\}}.$
(3.55)

  3. 3.

$\hat{q}^{(r)}(u_{i})=\mathrm{Beta}\left(\omega^{(r)}_{i},\nu^{(r)}_{i}\right)$
for all $i\in\\{1,\dots,L\\}$ where $\omega^{(r)}_{i}$ and $\nu^{(r)}_{i}$ are
defined as:

$\displaystyle\omega^{(r)}_{i}=\delta_{u}\left(\omega^{(r-1)}_{i}-1\right)+\delta_{z}\sum_{j\in\mathcal{V}}\tau_{ji}^{(r)}+1,$
$\displaystyle\nu^{(r)}_{i}=\delta_{u}\left(\nu^{(r-1)}_{i}-1\right)+\delta_{z}\sum_{j\in\mathcal{V}}\sum_{k=i+1}\tau_{jk}^{(r)}+1.$
(3.56)

A subtle problem arises with the updates in (3.53). In the case that the
truncation parameter, $L$, is larger than the true number of groups, the
update procedure will allocate some groups zero nodes. If
$\ell\in\\{1,\cdots,L\\}$ is one such group, then for all $i\in\mathcal{V}$
$\tau_{i\ell}^{(r)}=0$ for every $r$ after the algorithm has converged, up
until a change occurs. For $\delta_{\lambda}\in(0,1)$, it follows from the
rate update in (3.53) that $\lim_{r\to\infty}\beta_{\ell
m}^{(r)}=\lim_{r\to\infty}\beta_{m\ell}^{(r)}=0$ for all
$m\in\\{1,\cdots,L\\}$, and thus that the posterior mean of $\lambda_{\ell m}$
and $\lambda_{m\ell}$ diverges. This divergence causes problems in (3.55)
where this mean appears. To circumvent this problem, we must replace
$\delta_{\lambda}$ with the set $\\{\delta^{km}_{\lambda}\\}_{k,m=1}^{L}$.
That is, we introduce a specific forgetting factor for each group-to-group
rate. For each, $k,m\in\\{1,\dots,L\\}$, we initially set
$\delta^{km}_{L}\equiv\delta_{\lambda}$, but monitor the sum
$\sum_{(i,j)\in\mathcal{E}}\tau_{ik}^{(r-1)}\tau_{jm}^{(r-1)}$ with increasing
$r$. Then, for all $k,m\in\\{1,\dots,L\\}$ such that
$\sum_{(i,j)\in\mathcal{E}}\tau_{ik}^{(r-1)}\tau_{jm}^{(r-1)}<\epsilon$, we
set $\delta_{km}=1$, where $\epsilon$ is some threshold. This intervention
prevents the exponential decay of $\beta_{km}^{(r)}$ to 0 as $r\to\infty$.
Once the threshold is exceeded, the relevant BFFs are returned to
$\delta_{\lambda}$. Experimentation found that $0.1$ is a good choice for
$\epsilon$.

Figure 8: A directed acyclic graph of the model given by (3.43)-(3.46).

##  4 Simulation studies

We evaluate the online changepoint detection algorithm of Section 3 using
simulated data. Four simulation studies are conducted. Unless otherwise
stated, we consider a network with $N=500$ nodes and $K=2$ latent groups,
initialised with $\Pi=(0.6,0.4)$. For convenience, we define the intra-inter-
group rate matrix

$\lambda_{0}=\begin{pmatrix}2&1\\\ 0.3&8\end{pmatrix}.$ (4.1)

The update interval is set to $\Delta=0.1$ time units throughout, and every
experiment is repeated 50 times, with the results averaged. We cycle 3 times
over the CAVI and fixed point equations. All hyperparameters are initialised
as 1, except for $\gamma$ and $\xi$, which are initialised uniformly on
$[0.95,1.05]^{K}$ and $[0.95,1.05]^{K^{\prime}}$, respectively. For all
$i,j\in\mathcal{V}$, $\sigma_{ij}$ is initialised as $1/2$, and $\tau_{i}$ as
$\bm{1}_{K}/K$, when $K$ is known, and as $\bm{1}_{L}/L$ when we infer the
number of groups at a truncation level of $L$. Unless otherwise stated,
$L_{JS}$ is set to be 2 and $L_{KL},B_{1}$ and $B_{2}$ to 10. In Section 4.1,
we examine membership and rate recovery in the case that $\bm{A}$ and $K$ are
known when a varying proportion of nodes swap between the two groups. In
Section 4.2, membership recovery is investigated when $\bm{A}$ is known, but
$K$ changes throughout the observation window. In Section 4.3, we consider two
the effect of a decreasing lag between two instantaneous changes to $\lambda$
when both $\bm{A}$ and $K$ are known. Finally, in Section 4.4, the effect of
increasing sparsity in case of known $K$ but unknown $\bm{A}$ is examined.

###  4.1 Latent group membership recovery

The first experiment takes a fully-connected network with rate matrix
$\lambda_{0}$ as given in (4.1). At time $t=3$, $P\%$ of nodes from group 1
swap to group 2, with $P\in\\{1,10,25,50,75,95\\}$. To evaluate latent
membership recovery, the inferred group memberships at time $t$ are compared
to the true memberships using the adjusted rand index (ARI, Hubert and Arabie,
1985). The ARI takes values between $-1$ and $1$, with a value of $1$
indicating perfect agreement (up to label switching), $0$ random agreement,
and $-1$ complete disagreement. By computing the ARI at each update, rather
than taking an average over $[0,t]$, we can examine the smoothness of
recovery, and the reaction of the algorithm to changepoints.

(a) Mean ARI of the repetitions against update time for a varying proportion
of group 1 nodes changing to group 2.

(b) Mean over the repetitions of the proportion of nodes in group 1 with time.
The results obtained when using no BFF (green line) and a BFF of 0.1 (pink
line) are shown together on each panel. The true proportion of group 1 either
side of the changepoint is plotted as black, horizontal lines..

Figure 9: Detection of group membership changes for a varying percentage of
nodes swapping from group 1 to group 2 at $t=3$. The panel titles give the
percentage of nodes that swap from group 1 to 2.

Figure 9(a) shows the inference procedure is stable across all values of $P$,
where the plots begin from $B_{1}$ steps. The ARI is steady before and after
the change except for the step immediately after the change. Figure 9(b) shows
that the algorithm quickly converges to the true proportions, both with and
without a BFF.

###  4.2 Number of groups recovery

We consider a fully connected network with two group membership changes: all
nodes merging into group 1 at $t=2.5$, followed by the creation of a new group
at $t=3.5$. Specifically, at $t=3.5$, $P\%$ of the nodes in group 1 remain,
while the remainder create group 2. We again consider
$P\in\\{1,10,25,50,75,95\\}$, with $\lambda_{0}$ as the rate matrix, and ARI
as our performance metric. We set $L_{JS}=0$ here, which corresponds to taking
the argmax of each $\tau_{i}$.

Figure 11 shows that the inference correctly groups the nodes in all cases
pre-merger and between the merger and creation, both with and without a BFF.
Using a BFF, the algorithm is seen to maintain a high ARI after the creation
of a new group, in all cases except when $P=75$. The application without the
BFF also fails in this case, but additionally it fails for $P=50$, and is
slower to converge to the new groups when $P=25$ and 10.

For cases where 50% of nodes or less remain in group 1, the algorithm
converges to the correct proportions, but with label-switching. This is
observed in Figure 11, where the rates are seen to swap in the 10% panel for
the case of a BFF.

Figure 10: Mean ARI against update time for the merger of group 2 into 1 at
$t=2.5$, and the creation of group 2 at $t=3.5$, with the panel titles giving
the percentage of nodes remaining in group 1 after $t=3.5$. The black, dashed
vertical lines mark the changepoints. At $t=2.5$, all nodes in group 2 change
to group 1, and at $t=3.5$, $P\%$ of group 1 nodes change to group 2,
$P\in\\{1,10,25,50,75,95\\}$.

Figure 11: Posterior means and 95% confidence interval for components
$\lambda$ with update time. The black, dashed vertical lines mark the
changepoints. At $t=2.5$, all nodes in group 2 change to group 1, and at
$t=3.5$, $P\%$ of group 1 nodes change to group 2, $P\in\\{10,75,95\\}$.

###  4.3 Detection of changes to the matrix of rates

On a fully connected network, we examine sequential changes to the rate matrix
with decreasing time between the changes. The rate matrix maps as

$\lambda=\begin{pmatrix}2&1\\\
0.3&8\end{pmatrix}\quad\mapsto\quad\lambda^{\prime}=\begin{pmatrix}5&1\\\
0.3&8\end{pmatrix}\quad\mapsto\quad\lambda^{\prime\prime}=\begin{pmatrix}3&1\\\
0.3&8\end{pmatrix},$

with the first change at $t=3$ and the second at $t=3+0.1M$. For each run, $M$
took a value in $\\{1,2,3,4,5,10\\}$. Here we do not reset the stream after a
flagged change to the latent rates as we control their relative magnitude and
wish to examine the performance of the algorithm with small latency between
changes.

We evaluate the detection of latent changes using the proportion of
changepoints correctly detected (CCD), and the proportion of detections that
are not false (DNF), as in Bodenham and Adams (2017). Note that we aggregate
across all group-to-group rates. Specifically, if there are $C$ changepoints,
and we detect $D$, $T$ of which are correct, analogous to recall and
precision, the authors define:

  1. 1.

$\mathrm{CCD}=T/C$, the proportion of changepoints correctly detected,

  2. 2.

$\mathrm{DNF}=T/D$, the proportion of detections that are not false.

As noted by Bodenham and Adams (2017), CCD and DNF are preferred to their more
intuitive counterparts (proportion of missed changepoints and proportion of
false detections, respectively) as for the CCD and DNF, values closer to 1
indicate better performance than those closer to 0. Furthermore, these metrics
are preferable over average run length metrics ($\mathrm{ARL}0$ and
$\mathrm{ARL}1$) as we want to capture the number of changes detected and
missed.

Suppose $t^{\ast}_{n}$ and $t^{\ast}_{n+1}$ are consecutive true changepoints
for $\lambda_{km}$, and let the most recently flagged changepoint by the
algorithm for $\lambda_{km}$ be $t^{\prime}_{\ell}<t^{\ast}_{n}$. If the next
flag for $\lambda_{km}$ occurs at $t^{\prime}_{\ell+1}>t^{\ast}_{n+1}$, that
is $t^{\ast}_{n}$ is missed, then $t^{\prime}_{m+1}$ is classified as a
correct detection of $t^{\ast}_{n+1}$. This is common practice in the
literature. This approach is adopted by Bodenham and Adams (2017), whereas
other authors implement a softmax rule for classification (see, for example,
Alanqary et al., 2021; Yamanishi and Takeuchi, 2002; Bodenham and Adams,
2017), adjusting for when multiple changes are flagged in the same window.

Figure 12: Boxplots of CCD and DNF over 50 runs for no BFF and a BFF of
$\delta=0.1$. Each simulation has one change at $t=3$ and another a
$t=3+0.1M$, where $M$ is on the horizontal axis.

In Figure 12, we see the CCD and DNF for each run. The CCD is consistently
high for both a BFF and no BFF. This demonstrates that the algorithm
consistently flags true changes with and without a BFF. The role of the BFF is
shown in the DNF, where we see BFF yielding higher values. Without a BFF,
multiple changes are flagged that did not occur.

###  4.4 Effect of network sparsity

We examine the effect of decreasing network sparsity on changepoint detection.
We simulate a network with $K^{\prime}=1$ latent connection groups, letting
$\rho$, the group connection probability, vary over
$\\{0.01,0.025,0.05,0.1,0.25,0.5\\}$. For the point process groups, we retain
the same $\Pi$ as before, and again, we set $\lambda=\lambda_{0}$. We simulate
on $[0,25]$, with $\Delta=0.1$, and at $t=10$, $25\%$ of nodes from group 1
swap to group 2. The inference procedure is run on the simulated network
twice, once where the network is assumed fully-connected, and a second time
when the adjacency matrix is inferred. Both procedures are run with
$\delta=0.1$.

Figure 14 demonstrates the effect of incorrectly assuming that the graph is
fully connected: the rates are significantly underestimated. On the other
hand, when $\bm{A}$ is inferred, the posterior means are much closer to the
true values, even in the case of $1\%$ density. Furthermore, in Figure 14 the
mean ARI of the membership recovery is seen to be higher for the inferred
adjacency matrix.

Figure 13: Posterior means and 95% confidence interval for components
$\lambda$ with update time. The black, dashed vertical line marks the
changepoint, at which $25\%$ of nodes in group 1 swap to group 2.

Figure 14: Mean ARI of the repetitions against update time. The black, dashed
vertical line marks the changepoint, at which $25\%$ of nodes in group 1 swap
to group 2.

##  5 Application to the Santander Cycles bike-sharing network

The proposed online changepoint algorithms were tested on Transport for London
(TfL) data from the London Santander Cycles bike-sharing network, which is
publicly available online (https://cycling.data.tfl.gov.uk/, powered by TfL
Open Data). Each datum corresponds to a bike hire, and contains the start and
end times of the journey, the IDs of the source and destination stations, the
journey duration, a bike ID number, and an unique identifier for the journey.
Considering the start and end stations as source and destination nodes, and
the timestamp of the end of the journey as an arrival time to the directed
edge from source to destination, the data forms a network point process. In
this study, the data is aggregated into weekly counts to smooth the
intensities of the point processes and weekly periodicities. We select a
subset of the data from 2nd January 2019 until 15th July 2020 as this window
contains significant COVID-19 related national events that can be used to
check the performance of our algorithm. In this time period, $N=791$ unique
nodes are observed within $T=80$ weekly time windows, with updates every
$\Delta=1$ week time steps.

The online VB algorithm for the dynamic BHPP is run for the separate cases of
an unknown graph structure, and an unknown number of groups. The number of
groups in the case of an unknown adjacency matrix was set to match the number
inferred by the implementation for an unknown number of groups, which was
$K=6$. The adjacency matrix in the case of unknown $K$ was set to correspond
to a fully connected graph: $\bm{A}=\bm{1}_{791\times 791}$. We consider the
task of detecting changes to the latent group structure of the bike sharing
network. We run both algorithms with $L_{JS}=1.55$, $B_{1}=25$, $B_{2}=10$,
and $\kappa=2$, so that changes can be detected only after
$B_{1}+B_{2}+\kappa=37$ weeks of observations. We initialise the
hyperparameters as described in Section 4.

Figure 15: The number of nodes that change memberships at each update point
from initialisation. The green regions correspond to the Christmas and New
Year period of 2019, the introduction of the first UK COVID-19 lockdown, and
the phased easing of these restrictions, respectively.

Figure 15 shows that the algorithms flag multiple changes at each update, but
that there are three regions where the number of flagged changes peaks. From
left to right, these green regions correspond to the Christmas and New Year
period of 2019, the introduction of the first UK COVID-19 lockdown on
19/03/2020, and the subsequent phased easing of restrictions from 01/06/2020.
The algorithm reacts to these events, which are likely to cause changes to the
network.

Figure 16: Station locations with those stations that change group at the
onset of the UK COVID-19 lockdown coloured dark red.

Figure 16 shows that the flagged changes corresponding to the onset of the
lockdown are very concentrated around central London, which makes sense as
“work from home” orders will have affected the use of these commuting bike
stations. Additionally, stations around the Westfield Shopping Centre in
Shepherd’s Bush and the Canary Wharf financial district are flagged,
representing an expected change due to the government restrictions.

Furthermore, Figure 17 displays an example clustering from the dynamic BHPP
model with unknown number of groups at update point 50. There is some spatial
clustering to the nodes, but each cluster contains mainly nodes with similar
activity patterns. For example, the dark and light blue clusters represent
popular nodes in central London and the Canary Wharf financial district, which
are mainly used by commuters into these areas. Similarly, the light and dark
green clusters are nodes in West and East London (with the dark green cluster
mostly covering West London, and the light green East London and the areas of
Battersea and Wandsworth (south of the river Thames). The pink cluster is
around the boundary between the blue and green regions. A notable exception is
represented by the red cluster, representing the most popular stations within
the network, mainly used for leisure around Hyde Park and the Queen Elizabeth
Olympic Park in Stratford. In addition, the red cluster also contains two
additional stations in Battersea Park and Ravenscourt Park (west of
Hammersmith).

Figure 17: Station locations coloured by assigned group at update step 50 of
the dynamic BHPP with an unknown number of groups and
$\bm{A}=\bm{1}_{791\times 791}$.

##  6 Conclusion and discussion

We have presented a novel online Bayesian inference framework for detecting
changes to the latent structure of a block-homogeneous Poisson process in
which data arrives as batches in an event stream. Our methodology is scalable,
and leverages a Bayesian forgetting factor framework to flag changes to the
latent community structure and edge-process rates. The framework is extended
to the cases where the adjacency matrix or the number of latent groups are
unknown a priori. When tested on both real and simulated data, our methodology
is seen to detect latent structure accurately and with minimal latency.

There are numerous ways in which this work could be extended. In particular,
the frameworks for an unknown graph structure and unknown number of latent
groups can readily be integrated to handle the case where neither is known a
priori. It would also be of interest to incorporate seasonality into the
model, perhaps within an online framework with longer memory.

A further challenge would be to adapt the framework to allow for nodes to
enter or leave the network during observation. Similarly, the adjacency matrix
is assumed static in our methodology, and so an extension would allow for this
to be dynamic, although this would likely cause identifiability issues in the
case where there are also changes to the latent rates.

## Acknowledgements

Joshua Corneck acknowledges funding from the Engineering and Physical Sciences
Research Council (EPSRC), grant number EP/S023151/1. Ed Cohen acknowledges
funding from the EPSRC NeST Programme, grant number EP/X002195/1. Francesco
Sanna Passino acknowledges funding from the EPSRC, grant number EP/Y002113/1.

## Code

Python code to implement the methodologies proposed in this article and
reproduce the results is available in the Github repository
joshcorneck/dynamicBHPP.

## References

  * Airoldi et al. (2013) Airoldi, E. M., Costa, T. B., and Chan, S. H. (2013) Stochastic blockmodel approximation of a graphon: Theory and consistent estimation.  In Advances in Neural Information Processing Systems (eds. C. Burges, L. Bottou, M. Welling, Z. Ghahramani and K. Weinberger), vol. 26. Curran Associates, Inc.
  * Alanqary et al. (2021) Alanqary, A., Alomar, A. O., and Shah, D. (2021) Change point detection via multivariate singular spectrum analysis.  In Advances in Neural Information Processing Systems (eds. A. Beygelzimer, Y. Dauphin, P. Liang and J. W. Vaughan).
  * Amini et al. (2013) Amini, A. A., Chen, A., Bickel, P. J., and Levina, E. (2013) Pseudo-likelihood methods for community detection in large sparse networks.  The Annals of Statistics, 41, 2097–2122.
  * Bifet et al. (2018) Bifet, A., Gavaldà, R., Holmes, G., and Pfahringer, B. (2018) Machine Learning for Data Streams: with Practical Examples in MOA.  Adaptive Computation and Machine Learning series. MIT Press.
  * Bishop (2006) Bishop, C. M. (2006) Pattern Recognition and Machine Learning (Information Science and Statistics).  Berlin, Heidelberg: Springer-Verlag.
  * Blei and Jordan (2006) Blei, D. M. and Jordan, M. I. (2006) Variational inference for Dirichlet process mixtures.  Bayesian Analysis, 1, 121–143.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017) Variational inference: A review for statisticians.  Journal of the American Statistical Association, 112, 859–877.
  * Bodenham and Adams (2017) Bodenham, D. A. and Adams, N. M. (2017) Continuous monitoring for changepoints in data streams using adaptive estimation.  Statistics and Computing, 27, 1257–1270.
  * Fang et al. (2024) Fang, G., Ward, O. G., and Zheng, T. (2024) Online estimation and community detection of network point processes for event streams.  Statistics and Computing, 34.
  * Ferguson (1973) Ferguson, T. S. (1973) A Bayesian analysis of some nonparametric problems.  The Annals of Statistics, 1, 209–230.
  * Gao et al. (2018) Gao, C., Ma, Z., Zhang, A. Y., and Zhou, H. H. (2018) Community detection in degree-corrected block models.  The Annals of Statistics, 46, 2153–2185.
  * Hallgren et al. (2023) Hallgren, K. L., Heard, N. A., and Turcotte, M. J. M. (2023) Changepoint detection on a graph of time series.  Bayesian Analysis, 1–28.
  * Haykin (2002) Haykin, S. (2002) Adaptive filter theory.  Upper Saddle River, NJ: Prentice Hall, 4th edn.
  * Heard et al. (2010) Heard, N. A., Weston, D. J., Platanioti, K., and Hand, D. J. (2010) Bayesian anomaly detection methods for social networks.  The Annals of Applied Statistics, 4, 645–662.
  * Holland et al. (1983) Holland, P. W., Laskey, K. B., and Leinhardt, S. (1983) Stochastic blockmodels: First steps.  Social Networks, 5, 109–137.
  * Holme (2015) Holme, P. (2015) Modern temporal network theory: a colloquium.  The European Physical Journal B, 88, 234.
  * Hubert and Arabie (1985) Hubert, L. J. and Arabie, P. (1985) Comparing partitions.  Journal of Classification, 2, 193–218.
  * Jasra et al. (2005) Jasra, A., Holmes, C. C., and Stephens, D. A. (2005) Markov Chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling.  Statistical Science, 20, 50–67.
  * Karrer and Newman (2011) Karrer, B. and Newman, M. E. J. (2011) Stochastic blockmodels and community structure in networks.  Physical Review E, 83, 016107.
  * Klimt and Yang (2004) Klimt, B. and Yang, Y. (2004) The Enron corpus: A new dataset for email classification research.  In Machine Learning: ECML 2004 (eds. J.-F. Boulicaut, F. Esposito, F. Giannotti and D. Pedreschi), 217–226. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Lee and Wilkinson (2019) Lee, C. and Wilkinson, D. J. (2019) A review of stochastic block models and extensions for graph clustering.  Applied Network Science, 4, 1–50.
  * Leys et al. (2013) Leys, C., Ley, C., Klein, O., Bernard, P., and Licata, L. (2013) Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median.  Journal of Experimental Social Psychology, 49, 764–766.
  * Lin (1991) Lin, J. (1991) Divergence measures based on the Shannon entropy.  IEEE Transactions on Information Theory, 37, 145–151.
  * Mariadassou et al. (2010) Mariadassou, M., Robin, S., and Vacher, C. (2010) Uncovering latent structure in valued graphs: A variational approach.  The Annals of Applied Statistics, 4.
  * Mariadassou and Tabouy (2020) Mariadassou, M. and Tabouy, T. (2020) Consistency and asymptotic normality of stochastic block models estimators from sampled data.  The Electronic Journal of Statistics, 14, 3672–3704.
  * Matias and Miele (2017) Matias, C. and Miele, V. (2017) Statistical clustering of temporal networks through a dynamic stochastic block model.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 79, 1119–1141.
  * Matias et al. (2018) Matias, C., Rebafka, T., and Villers, F. (2018) A semiparametric extension of the stochastic block model for longitudinal networks.  Biometrika, 105, 665–680.
  * Paranjape et al. (2017) Paranjape, A., Benson, A. R., and Leskovec, J. (2017) Motifs in temporal networks.  In Proceedings of the Tenth ACM International Conference on Web Search and Data Mining, WSDM ’17, 601–610. New York, NY, USA: Association for Computing Machinery.
  * Perry and Wolfe (2013) Perry, P. O. and Wolfe, P. J. (2013) Point process modelling for directed interaction networks.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 75, 821–849.
  * Pitman (2002) Pitman, J. (2002) Combinatorial stochastic processes (Technical Report 621). Department of Statistics.  University of California at Berkeley, Berkeley, CA.
  * Ray and Szabó (2022) Ray, K. and Szabó, B. (2022) Variational Bayes for high-dimensional linear regression with sparse priors.  Journal of the American Statistical Association, 117, 1270–1281.
  * Sanna Passino and Heard (2023) Sanna Passino, F. and Heard, N. A. (2023) Mutually exciting point process graphs for modeling dynamic networks.  Journal of Computational and Graphical Statistics, 32, 116–130.
  * Sethuraman (1994) Sethuraman, J. (1994) A constructive definition of Dirichlet priors.  Statistica Sinica, 639–650.
  * Shlomovich et al. (2022) Shlomovich, L., Cohen, E. A. K., and Adams, N. (2022) A parameter estimation method for multivariate binned Hawkes processes.  Statistics and Computing, 32, 98.
  * van Dyk and Park (2008) van Dyk, D. A. and Park, T. (2008) Partially collapsed Gibbs samplers.  Journal of the American Statistical Association, 103, 790–796.
  * Wang and Blei (2019) Wang, Y. and Blei, D. M. (2019) Frequentist consistency of variational Bayes.  Journal of the American Statistical Association, 114, 1147–1161.
  * Xu and Hero (2014) Xu, K. and Hero, A. (2014) Dynamic stochastic blockmodels for time-evolving social networks.  IEEE Journal of Selected Topics in Signal Processing, 8, 552–562.
  * Yamanishi and Takeuchi (2002) Yamanishi, K. and Takeuchi, J. (2002) A unifying framework for detecting outliers and change points from non-stationary time series data.  In Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’02, 676–681. New York, NY, USA: Association for Computing Machinery.
  * Yang et al. (2011) Yang, T., Chi, Y., Zhu, S., Gong, Y., and Jin, R. (2011) Detecting communities and their evolutions in dynamic social networks - a Bayesian approach.  Machine Learning, 82, 157–189.
  * Zhang and Zhou (2020) Zhang, A. Y. and Zhou, H. H. (2020) Theoretical and computational guarantees of mean field variational inference for community detection.  The Annals of Statistics, 48, 2575–2598.

##  Appendix A Derivation of CAVI updates

Here we derive the CAVI approximating distributions for each parameter at the
$r$th update. In deriving these expressions, we assume no ordering for the
parameter updates, and will use a superscript of $(r-1)$ for all parameters
but for the one whose expression is then being derived.

###  A.1 Approximation of $q(z_{i})$

Ignoring any terms that do not contain $z_{i}$, we can derive the expectation
as follows:

$\displaystyle\mathbb{E}_{-z_{i}}$ $\displaystyle\\{\log
p(x,z,\lambda,\pi)\\}=\mathbb{E}_{-z_{i}}\left\\{\log
p(x|z,\lambda)+\delta_{z}\log p(z|\pi)\right\\}+\text{cst.}$ (A.1)
$\displaystyle=$ $\displaystyle\
\mathbb{E}_{-z_{i}}\Bigg{[}\sum_{j:(i,j)\in\mathcal{E}}\sum_{k,m\in\mathcal{K}}\tilde{z}_{ik}\tilde{z}_{jm}\bigg{(}x_{ij}(I_{r})\log(\lambda_{km})-\Delta\lambda_{km}\bigg{)}+$
(A.2)
$\displaystyle\qquad\qquad\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\sum_{k,m\in\mathcal{K}}\tilde{z}_{j^{\prime}k}\tilde{z}_{im}\bigg{(}x_{j^{\prime}i}(I_{r})\log(\lambda_{km})-\Delta\lambda_{km}\bigg{)}-$
(A.3)
$\displaystyle\qquad\qquad\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\bigg{(}x_{ii}(t)\log(\lambda_{kk})-\Delta\lambda_{kk}\bigg{)}+\delta_{z}\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\log\pi_{k}\Bigg{]}+\text{cst.}$
(A.4) $\displaystyle=$
$\displaystyle\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\Bigg{[}\delta_{z}\mathbb{E}_{\pi}\\{\log\pi_{k}\\}+x_{ii}(t)\mathbb{E}_{\lambda_{kk}}\\{\log(\lambda_{kk})\\}-\Delta\mathbb{E}_{\lambda_{kk}}\\{\lambda_{kk}\\}+$
(A.5)
$\displaystyle\sum_{m\in\mathcal{K}}\Bigg{\\{}\sum_{j:(i,j)\in\mathcal{E}}\mathbb{E}_{z_{j}}\\{\tilde{z}_{jm}\\}\bigg{(}x_{ij}(I_{r})\mathbb{E}_{\lambda_{km}}\\{\log(\lambda_{km})\\}-\Delta\mathbb{E}_{\lambda_{km}}\\{\lambda_{km}\\}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}+$
(A.6)
$\displaystyle\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\mathbb{E}_{\tilde{z}_{j^{\prime}}}\\{\tilde{z}_{j^{\prime}m}\\}\bigg{(}x_{j^{\prime}i}(I_{r})\mathbb{E}_{\lambda_{mk}}\\{\log(\lambda_{mk})\\}-\Delta\mathbb{E}_{\lambda_{mk}}\\{\lambda_{mk}\\}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j^{\prime})\bigg{)}\Bigg{\\}}\Bigg{]}+\text{cst.}$
(A.7)

Taking the exponential, it follows that the optimal choice of $\hat{q}(z_{i})$
under a CAVI approximation is
$\hat{q}^{(r)}(z_{i})=\mathrm{Categorical}\left(z_{i};\tau_{i}^{(r)}\right)$,
where

$\displaystyle\tau_{ik}^{(r)}$
$\displaystyle\propto\exp\Bigg{\\{}\delta_{z}\Bigg{[}\psi\left(\gamma^{(r-1)}_{k}\right)-\psi\left(\sum_{\ell}\gamma^{(r-1)}_{\ell}\right)\Bigg{]}+x_{ii}(t)\bigg{(}\psi\left(\alpha^{(r-1)}_{kk}\right)-\log\left(\beta^{(r-1)}_{kk}\right)\bigg{)}-\Delta\frac{\alpha^{(r-1)}_{kk}}{\beta^{(r-1)}_{kk}}+$
(A.8)
$\displaystyle\qquad\sum_{m\in\mathcal{K}}\Bigg{[}\sum_{j:(i,j)\in\mathcal{E}}\tau^{(r-1)}_{jm}\bigg{(}x_{ij}(I_{r})\left(\psi\left(\alpha^{(r-1)}_{km}\right)-\log\left(\beta^{(r-1)}_{km}\right)\right)-\Delta\frac{\alpha^{(r-1)}_{km}}{\beta^{(r-1)}_{km}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}\
+$ (A.9)
$\displaystyle\qquad\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\tau_{j^{\prime}m}^{(r-1)}\bigg{(}x_{j^{\prime}i}(I_{r})\left(\psi\left(\alpha^{(r-1)}_{mk}\right)-\log\left(\beta^{(r-1)}_{mk}\right)\right)-\Delta\frac{\alpha^{(r-1)}_{mk}}{\beta^{(r-1)}_{mk}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j^{\prime})\bigg{)}\Bigg{]}\Bigg{\\}},$
(A.10)

for $i\in\mathcal{V}$ and $k\in\mathcal{K}$.

###  A.2 Approximation of $q(\pi)$

Ignoring any terms that do not contain $\pi$, we compute the expectation as:

$\displaystyle\mathbb{E}_{-\pi}$ $\displaystyle\left\\{\log
p(x,z,\lambda,\pi)\right\\}=\mathbb{E}_{-\pi}\left\\{\delta_{z}\log
p(z|\pi)+\delta_{\pi}\log p(\pi)\right\\}+\text{cst.}$ (A.11)
$\displaystyle=\sum_{k\in\mathcal{K}}\log\pi_{k}\left[\sum_{i\in\mathcal{V}}\delta_{z}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}+\delta_{\pi}\left(\gamma_{k}^{(r-1)}-1\right)\right]+\text{cst.}$
(A.12)

Using the distributions derived for $\hat{q}^{(r-1)}(z_{i})$, we see that the
optional CAVI approximation is to take
$\hat{q}^{(r)}(\pi)=\mathrm{Dirichlet}\left(\pi;\gamma^{(r)}\right)$, where

$\gamma_{k}^{(r)}=\delta_{\pi}\left(\gamma_{k}^{(r-1)}-1\right)+\delta_{z}\sum_{i\in\mathcal{V}}\tau^{(r-1)}_{ik}+1,$
(A.13)

for $k\in\mathcal{K}$.

###  A.3 Approximation of $q(\lambda_{km})$

Ignoring any terms that do not contain $\lambda_{km}$, we derive the
expectation as follows:

$\displaystyle\mathbb{E}_{-\lambda_{km}}$ $\displaystyle\left\\{\log
p(x,z,\lambda,\pi)\right\\}=\mathbb{E}_{-\lambda_{km}}\left\\{\log
p(x|z,\lambda)+\delta_{\lambda}\log p(\lambda)\right\\}+\text{cst.}$ (A.14)
$\displaystyle=\
\sum_{(i,j)\in\mathcal{E}}\sum_{k,m\in\mathcal{K}}\bigg{[}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{j}}\\{\tilde{z}_{jm}\\}(x_{ij}\log\lambda_{km}-\lambda_{km}\Delta)\bigg{]}+$
(A.15)
$\displaystyle\quad\delta_{\lambda}\sum_{k,m\in\mathcal{K}}\left[\left(\alpha^{(r-1)}_{km}-1\right)\log\lambda_{km}-\beta^{(r-1)}_{km}\lambda_{km}\right]+\text{cst.}$
(A.16)
$\displaystyle=\sum_{k,m\in\mathcal{K}}\bigg{\\{}\bigg{(}\delta_{\lambda}\left(\alpha^{(r-1)}_{km}-1\right)+\sum_{(i,j)\in\mathcal{E}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{j}}\\{\tilde{z}_{ik}\\}x_{ij}(I_{r})\bigg{)}\log\lambda_{km}\
-$ (A.17)
$\displaystyle\quad\left(\delta_{\lambda}\beta^{(r-1)}_{km}+\Delta\sum_{(i,j)\in\mathcal{E}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{j}}\\{\tilde{z}_{ik}\\}\right)\lambda_{km}\bigg{\\}}+\text{cst.}$
(A.18)

Using the distributions derived for $\hat{q}^{(r-1)}(z_{i})$, it follows that
the optimal CAVI approximation is to take
$\hat{q}^{(r)}\left(\lambda_{km}\right)=\mathrm{Gamma}\left(\lambda_{km};\alpha^{(r)}_{km},\beta^{(r)}_{km}\right)$,
where

$\displaystyle\alpha^{(r)}_{km}$
$\displaystyle=\delta_{\lambda}\left(\alpha^{(r-1)}_{km}-1\right)+\sum_{(i,j)\in\mathcal{E}}\tau^{(r-1)}_{ik}\tau^{(r-1)}_{jm}x_{ij}(I_{r})+1,$
(A.19) $\displaystyle\beta^{(r)}_{km}$
$\displaystyle=\delta_{\lambda}\beta^{(r-1)}_{km}+\Delta\sum_{(i,j)\in\mathcal{E}}\tau^{(r-1)}_{ik}\tau^{(r-1)}_{jm},$
(A.20)

for $k,m\in\mathcal{K}$.

##  Appendix B Derivation of CAVI updates for an unknown adjacency matrix

Here we derive the CAVI approximating distributions for each parameter at the
$r$th update. In deriving these expressions, we assume no ordering for the
parameter updates, and will use a superscript of $(r-1)$ for all parameters
but for the one whose expression is then being derived.

###  B.1 Approximation of $q(z_{i})$

Ignoring any terms that do not contain $z_{i}$, we derive the expectation as:

$\displaystyle\mathbb{E}_{-z_{i}}$ $\displaystyle\left\\{\log
p(x,z,a,\lambda,\pi)\right\\}=\mathbb{E}_{-z_{i}}\left\\{\log
p(x|z,a,\lambda)+\delta_{z}\log p(z|\pi)\right\\}+\text{cst.}$ (B.1)
$\displaystyle=\mathbb{E}_{-z_{i}}\Bigg{\\{}\sum_{j\in\mathcal{V}}\sum_{k,m\in\mathcal{K}}\tilde{z}_{ik}\tilde{z}_{jm}\bigg{[}a_{ij}\left(x_{ij}(I_{r})\log\lambda_{km}-\Delta\lambda_{km}\right)\bigg{]}+$
(B.2)
$\displaystyle\qquad\qquad\sum_{j^{\prime}\in\mathcal{V}}\sum_{k,m\in\mathcal{K}}\tilde{z}_{j^{\prime}k}\tilde{z}_{im}\bigg{[}a_{j^{\prime}i}\left(x_{j^{\prime}i}(I_{r})\log\lambda_{km}-\Delta\lambda_{km}\right)\bigg{]}-$
(B.3)
$\displaystyle\qquad\qquad\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\bigg{[}a_{ii}\left(x_{ii}(I_{r})\log\lambda_{kk}-\Delta\lambda_{kk}\right)\bigg{]}+\delta_{z}\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\log\pi_{k}\Bigg{\\}}+\text{cst.}$
(B.4)
$\displaystyle=\mathbb{E}_{-z_{i}}\Bigg{\\{}\sum_{j\in\mathcal{V}}\sum_{k,m\in\mathcal{K}}\tilde{z}_{ik}\tilde{z}_{jm}\bigg{[}a_{ij}\left(x_{ij}(I_{r})\log\lambda_{km}-\Delta\lambda_{km}\right)+a_{ji}\left(x_{ji}(I_{r})\log\lambda_{mk}-\Delta\lambda_{mk}\right)\bigg{]}-$
(B.5)
$\displaystyle\qquad\qquad\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\bigg{[}a_{ii}\left(x_{ii}(I_{r})\log\lambda_{kk}-\Delta\lambda_{kk}\right)\bigg{]}+\delta_{z}\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\log\pi_{k}\Bigg{\\}}+\text{cst.}$
(B.6)
$\displaystyle=\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\Bigg{\\{}\sum_{j\in\mathcal{V}}\sum_{m\in\mathcal{K}}\tilde{z}_{jm}\bigg{[}x_{ij}(I_{r})\mathbb{E}_{a_{ij}}\\{a_{ij}\\}\mathbb{E}_{\lambda_{km}}\\{\log\lambda_{km}\\}+x_{ji}(I_{r})\mathbb{E}_{a_{ji}}\\{a_{ji}\\}\mathbb{E}_{\lambda_{mk}}\\{\log\lambda_{mk}\\}\
-$ (B.7)
$\displaystyle\qquad\qquad\Delta\bigg{(}\mathbb{E}_{a_{ij}}\\{a_{ij}\\}\mathbb{E}_{\lambda_{km}}\\{\lambda_{km}\\}+\mathbb{E}_{a_{ji}}\\{a_{ji}\\}\mathbb{E}_{\lambda_{mk}}\\{\lambda_{mk}\\}\bigg{)}\bigg{]}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}+$
(B.8)
$\displaystyle\qquad\qquad\mathbb{E}_{a_{ii}}\\{a_{ii}\\}\bigg{(}x_{ii}(I_{r})\mathbb{E}_{\lambda_{kk}}\\{\log\lambda_{kk}\\}-\Delta\mathbb{E}_{\lambda_{kk}}\\{\lambda_{kk}\\}\bigg{)}+\delta_{z}\mathbb{E}_{\pi_{k}}\left\\{\log\pi_{k}\right\\}\Bigg{\\}}+\text{cst.}$
(B.9)

Taking the exponential, and using the distributions derived for
$\hat{q}^{(r-1)}(\lambda_{km})$ and the assumed form for
$\hat{q}^{(r-1)}(a_{ij})$, it follows that the optimal CAVI approximation is
to take
$\hat{q}^{(r)}(z_{i})=\mathrm{Categorical}\left(z_{i};\tau^{(r)}_{i}\right)$,
where

$\displaystyle\tau^{(r)}_{ik}$
$\displaystyle\propto\exp\Bigg{\\{}\sum_{j\in\mathcal{V}}\sum_{m\in\mathcal{K}}\tau^{(r-1)}_{jm}\bigg{[}\sigma^{(r-1)}_{ij}x_{ij}(I_{r})\bigg{(}\psi\left(\alpha^{(r-1)}_{km}\right)-\log\left(\beta^{(r-1)}_{km}\right)\bigg{)}\
+$ (B.10)
$\displaystyle\quad\sigma^{(r-1)}_{ji}x_{ji}(I_{r})\bigg{(}\psi\left(\alpha^{(r-1)}_{mk}\right)-\log\left(\beta^{(r-1)}_{mk}\right)\bigg{)}-\Delta\bigg{(}\sigma^{(r-1)}_{ij}\frac{\alpha^{(r-1)}_{km}}{\beta^{(r-1)}_{km}}+\sigma^{(r-1)}_{ji}\frac{\alpha^{(r-1)}_{mk}}{\beta^{(r-1)}_{mk}}\bigg{)}\bigg{]}\
\times$ (B.11)
$\displaystyle\quad\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}+\sigma^{(r-1)}_{ii}\bigg{(}x_{ii}(I_{r})\bigg{(}\psi\left(\alpha^{(r-1)}_{kk}\right)-\log\left(\beta^{(r-1)}_{kk}\right)\bigg{)}-\Delta\frac{\alpha^{(r-1)}_{kk}}{\beta^{(r-1)}_{kk}}\bigg{)}\
+$ (B.12)
$\displaystyle\quad\delta_{z}\Bigg{[}\psi\left(\gamma^{(r-1)}_{k}\right)-\psi\left(\sum_{\ell}\gamma^{(r-1)}_{\ell}\right)\Bigg{]}\Bigg{\\}},$
(B.13)

for $i\in\mathcal{V}$ and $k\in\mathcal{K}$.

###  B.2 Approximation of $q(\lambda_{km})$

Ignoring any terms that do not contain $\lambda_{km}$, we can derive the
expectation as follows:

$\displaystyle\mathbb{E}_{-\lambda_{km}}$ $\displaystyle\\{\log
p(x,z,z^{\prime},a,\lambda,\rho,\pi,\mu)\\}=\mathbb{E}_{-\lambda_{km}}\bigg{\\{}\log
p(x|z,a,\lambda)+\delta_{\lambda}\log p(\lambda)\bigg{\\}}+\text{cst.}$ (B.14)
$\displaystyle=\sum_{(i,j)\in\mathcal{R}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{j}}\\{\tilde{z}_{jm}\\}\mathbb{E}_{a_{ij}}\\{a_{ij}\\}\bigg{(}x_{ij}(I_{r})\log\lambda_{km}-\Delta\lambda_{km}\bigg{)}+$
(B.15)
$\displaystyle\quad\delta_{\lambda}\bigg{(}\bigg{(}\alpha_{km}^{(r-1)}-1\bigg{)}\log\lambda_{km}-\beta_{km}^{(r-1)}\lambda_{km}\bigg{)}+\text{cst.}$
(B.16)
$\displaystyle=\bigg{(}\sum_{(i,j)\in\mathcal{R}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{i}}\\{\tilde{z}_{jm}\\}\mathbb{E}+{a_{ij}}\\{a_{ij}\\}x_{ij}(I_{r})+\delta_{\lambda}\bigg{(}\alpha_{km}^{(r-1)}-1\bigg{)}\bigg{)}\log\lambda_{km}\
-$ (B.17)
$\displaystyle\quad\bigg{(}\Delta\sum_{(i,j)\in\mathcal{R}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{j}}\\{\tilde{z}_{jm}\\}\mathbb{E}_{a_{ij}}\\{a_{ij}\\}+\delta_{\lambda}\beta_{km}^{(r-1)}\bigg{)}\lambda_{km}+\text{cst.}$
(B.18)

Using the distributions derived for $\hat{q}^{(r-1)}(z_{i})$ and assumed form
for $\hat{q}^{(r-1)}(a_{ij})$, we see that the optimal CAVI approximating
distribution is
$\hat{q}^{(r)}(\lambda_{km})=\mathrm{Gamma}\bigg{(}\lambda_{km};\alpha_{km}^{(r)},\beta_{km}^{(r)}\bigg{)}$
with

$\displaystyle\alpha_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\bigg{(}\alpha_{km}^{(r-1)}-1\bigg{)}+\sum_{(i,j)\in\mathcal{R}}\tau^{(r-1)}_{ik}\tau^{(r-1)}_{jm}\sigma_{ij}^{(r-1)}x_{ij}(I_{r})+1,$
(B.19) $\displaystyle\beta_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\beta_{km}^{(r-1)}+\Delta\sum_{(i,j)\in\mathcal{R}}\tau^{(r-1)}_{ij}\tau^{(r-1)}_{jm}\sigma^{(r-1)}_{ij},$
(B.20)

for $k,m\in\mathcal{K}$.

###  B.3 Approximation of $q(\rho_{k^{\prime}m^{\prime}})$

Ignoring any terms that do not contain $\rho_{k^{\prime}m^{\prime}}$, we
derive the expectation as:

$\displaystyle\mathbb{E}_{-\rho_{k^{\prime}m^{\prime}}}$ $\displaystyle\\{\log
p(x,z,z^{\prime},a,\lambda,\rho,\pi,\mu)\\}=\mathbb{E}_{-\rho_{k^{\prime}m^{\prime}}}\\{\log
p(a|z,\rho)+\delta_{\rho}\log p(\rho)\\}+\text{cst.}$ (B.21)
$\displaystyle=\sum_{(i,j)\in\mathcal{R}}\mathbb{E}_{z_{i}^{\prime}}\\{\tilde{z}_{ik^{\prime}}^{\prime}\\}\mathbb{E}_{z_{j}^{\prime}}\\{\tilde{z}_{jm^{\prime}}^{\prime}\\}\bigg{[}\mathbb{E}_{a_{ij}}\\{a_{ij}\\}\log\rho_{k^{\prime}m^{\prime}}+\mathbb{E}_{a_{ij}}\\{1-a_{ij}\\}\log(1-\rho_{k^{\prime}m^{\prime}})\bigg{]}\
+$ (B.22)
$\displaystyle\quad\delta_{\rho}\bigg{(}\eta_{k^{\prime}m^{\prime}}^{(r-1)}-1\bigg{)}\log\rho_{k^{\prime}m^{\prime}}+\delta_{\rho}\bigg{(}\zeta_{k^{\prime}m^{\prime}}^{(r-1)}-1\bigg{)}\log(1-\rho_{k^{\prime}m^{\prime}})+\text{cst.}$
(B.23)
$\displaystyle=\bigg{[}\sum_{(i,j)\in\mathcal{R}}\mathbb{E}_{z_{i}^{\prime}}\\{\tilde{z}_{ik^{\prime}}^{\prime}\\}\mathbb{E}_{z_{j}^{\prime}}\\{\tilde{z}_{jm^{\prime}}^{\prime}\\}\mathbb{E}_{a_{ij}}\\{a_{ij}\\}+\delta_{\rho}\bigg{(}\eta_{k^{\prime}m^{\prime}}^{(r-1)}-1\bigg{)}\bigg{]}\log\rho_{k^{\prime}m^{\prime}}\
+$ (B.24)
$\displaystyle\quad\bigg{[}\sum_{(i,j)\in\mathcal{R}}\mathbb{E}_{z_{i}^{\prime}}\\{\tilde{z}_{ik^{\prime}}^{\prime}\\}\mathbb{E}_{z_{j}^{\prime}}\\{\tilde{z}_{jm^{\prime}}^{\prime}\\}\mathbb{E}_{a_{ij}}\\{1-a_{ij}\\}+\delta_{\rho}\bigg{(}\zeta_{k^{\prime}m^{\prime}}^{(r-1)}-1\bigg{)}\bigg{]}\log(1-\rho_{k^{\prime}m^{\prime}})+\text{cst.}$
(B.25)

Taking the exponential and computing expectations with respect to the derived
distributions for $\hat{q}^{(r-1)}(z_{i}^{\prime})$ and the assumed form for
$\hat{q}^{(r-1)}(a_{ij})$, it follows that the optimal CAVI approximation
takes
$\hat{q}^{(r)}(\rho_{k^{\prime}m^{\prime}})=\mathrm{Beta}\bigg{(}\rho_{k^{\prime}m^{\prime}};\eta^{(r)}_{k^{\prime}m^{\prime}},\zeta_{k^{\prime}m^{\prime}}^{(r)}\bigg{)}$,
with

$\displaystyle\eta^{(r)}_{k^{\prime}m^{\prime}}$
$\displaystyle=\delta_{\rho}\bigg{(}\eta_{k^{\prime}m^{\prime}}^{(r-1)}-1\bigg{)}+\sum_{(i,j)\in\mathcal{R}}\nu^{(r-1)}_{ik^{\prime}}\nu^{(r-1)}_{jm^{\prime}}\sigma_{ij}+1,$
(B.26) $\displaystyle\zeta^{(r)}_{k^{\prime}m^{\prime}}$
$\displaystyle=\delta_{\rho}\bigg{(}\zeta_{k^{\prime}m^{\prime}}^{(r-1)}-1\bigg{)}+\sum_{(i,j)\in\mathcal{R}}\nu^{(r-1)}_{ik^{\prime}}\nu^{(r-1)}_{jm^{\prime}}\left(1-\sigma^{(r-1)}_{ij}\right)+1,$
(B.27)

for $k^{\prime},m^{\prime}\in\mathcal{K}^{\prime}$.

###  B.4 Approximation of $q(\pi)$

Ignoring any terms that do not contain $\pi$, we derive the expectation as:

$\displaystyle\mathbb{E}_{-\pi}$ $\displaystyle\\{\log
p(x,z,z^{\prime},a,\lambda,\rho,\pi,\mu)\\}=\mathbb{E}_{-\pi}\\{\delta_{z}\log
p(z|\pi)+\delta_{\pi}\log p(\pi)\\}$ (B.28)
$\displaystyle=\mathbb{E}_{-\pi}\bigg{\\{}\delta_{z}\sum_{i\in\mathcal{V}}\sum_{k\in\mathcal{K}}\tilde{z}_{ik}\log\pi_{k}+\delta_{\pi}\sum_{k\in\mathcal{K}}\bigg{(}\gamma_{k}^{(r-1)}-1\bigg{)}\log\pi_{k}\bigg{\\}}+\text{cst.}$
(B.29)
$\displaystyle=\sum_{k\in\mathcal{K}}\bigg{\\{}\sum_{i\in\mathcal{V}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}+\delta_{\pi}\bigg{(}\gamma_{k}^{(r-1)}-1\bigg{)}\bigg{\\}}\log\pi_{k}+\text{cst.}$
(B.30)

Using the derived forms for $\hat{q}^{(r-1)}(z_{i})$, we see that the optimal
CAVI choice is
$\hat{q}^{(r)}(\pi)=\mathrm{Dirichlet}\left(\pi;\gamma^{(r)}\right)$, where

$\gamma^{(r)}_{k}=\delta_{\pi}\bigg{(}\gamma^{(r-1)}_{k}-1\bigg{)}+\delta_{z}\sum_{i\in\mathcal{V}}\tau^{(r-1)}_{ik}+1,$

for $k\in\mathcal{K}$.

###  B.5 Approximation of $q(\mu)$

Ignoring any terms that do not contain $\mu$, we compute the necessary
expectation as:

$\displaystyle\mathbb{E}_{-\mu}$ $\displaystyle\\{\log
p(x,z,z^{\prime},a,\lambda,\rho,\pi,\mu)\\}=\mathbb{E}_{-\mu}\\{\log
p(z^{\prime}|\mu)+\delta_{\mu}\log p(\mu)\\}$ (B.31)
$\displaystyle=\mathbb{E}_{-\mu}\bigg{\\{}\sum_{i\in\mathcal{V}}\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{ik^{\prime}}^{\prime}\log\mu_{k^{\prime}}+\delta_{\mu}\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\bigg{(}\xi_{k^{\prime}}^{(r-1)}-1\bigg{)}\log\mu_{k^{\prime}}\bigg{\\}}+\text{cst.}$
(B.32)
$\displaystyle=\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\bigg{\\{}\sum_{i\in\mathcal{V}}\mathbb{E}_{z_{i}^{\prime}}\\{\tilde{z}_{ik^{\prime}}^{\prime}\\}+\delta_{\mu}\bigg{(}\xi_{k^{\prime}}^{(r-1)}-1\bigg{)}\bigg{\\}}\log\mu_{k^{\prime}}+\text{cst.}$
(B.33)

Taking the exponential, and using the derived form for
$\hat{q}^{(r-1)}(z_{i}^{\prime})$, we see that the optimal CAVI distribution
is $\hat{q}^{(r)}(\mu)=\mathrm{Dirichlet}\left(\mu;\xi^{(r)}\right)$, where

$\xi^{(r)}_{k^{\prime}}=\delta_{\mu}\bigg{(}\mu^{(r-1)}_{k^{\prime}}-1\bigg{)}+\sum_{i\in\mathcal{V}}\nu^{(r-1)}_{ik^{\prime}}+1,$
(B.34)

for $k^{\prime}\in\mathcal{K}$.

###  B.6 Approximation of $q(z_{i}^{\prime})$

Ignoring terms that do not contain $z_{i}^{\prime}$, we derive the expectation
as:

$\displaystyle\mathbb{E}_{-z_{i}^{\prime}}\left\\{\log
p(x,z,z^{\prime},a,\lambda,\rho,\pi,\mu)\right\\}=\mathbb{E}_{-z_{i}^{\prime}}\left\\{\log
p(a|z^{\prime},\rho)+\log p(z^{\prime}|\mu)\right\\}+\text{cst.}$ (B.35)
$\displaystyle=\mathbb{E}_{-z_{i}^{\prime}}\Bigg{\\{}\sum_{i^{\prime},j^{\prime}\in\mathcal{V}}\sum_{k^{\prime},m^{\prime}\in\mathcal{K}^{\prime}}z_{ik^{\prime}}^{\prime}z_{jm^{\prime}}^{\prime}\bigg{[}a_{i^{\prime}j^{\prime}}\log\rho_{k^{\prime}m^{\prime}}+(1-a_{i^{\prime}j^{\prime}})\log(1-\rho_{k^{\prime}m^{\prime}})\bigg{]}+\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{i^{\prime}k^{\prime}}^{\prime}\log\mu_{k^{\prime}}\Bigg{\\}}+\text{cst.}$
(B.36)
$\displaystyle=\mathbb{E}_{-z_{i}^{\prime}}\Bigg{\\{}\sum_{j\in\mathcal{V}}\sum_{k^{\prime},m^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{jk^{\prime}}^{\prime}\tilde{z}_{im^{\prime}}^{\prime}\bigg{[}a_{ij}\log\rho_{k^{\prime}m^{\prime}}+(1-a_{ij})\log(1-\rho_{k^{\prime}m^{\prime}})\bigg{]}+$
(B.37)
$\displaystyle\qquad\qquad\sum_{j^{\prime}\in\mathcal{V}}\sum_{k^{\prime},m^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{ik^{\prime}}^{\prime}\tilde{z}_{j^{\prime}m^{\prime}}^{\prime}\bigg{[}a_{j^{\prime}i}\log\rho_{k^{\prime}m^{\prime}}+(1-a_{j^{\prime}i})\log(1-\rho_{k^{\prime}m^{\prime}})\bigg{]}-$
(B.38)
$\displaystyle\qquad\qquad\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{ik^{\prime}}^{\prime}\bigg{[}a_{ii}\log\rho_{k^{\prime}k^{\prime}}+(1-a_{ii})\log(1-\rho_{k^{\prime}k^{\prime}})\bigg{]}+\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{ik^{\prime}}^{\prime}\log\mu_{k^{\prime}}\Bigg{\\}}+\text{cst.}$
(B.39)
$\displaystyle=\sum_{k^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{ik^{\prime}}^{\prime}\Bigg{\\{}\sum_{j\in\mathcal{V}}\sum_{m^{\prime}\in\mathcal{K}^{\prime}}\tilde{z}_{jm^{\prime}}^{\prime}\bigg{[}\mathbb{E}_{a_{ij}}\\{a_{ij}\\}\mathbb{E}_{\rho_{k^{\prime}m^{\prime}}}\\{\log\rho_{k^{\prime}m^{\prime}}\\}+\mathbb{E}_{a_{ij}}\\{1-a_{ij}\\}\mathbb{E}_{\rho_{k^{\prime}m^{\prime}}}\\{\log(1-\rho_{k^{\prime}m^{\prime}})\\}\
+$ (B.40)
$\displaystyle\qquad\qquad\mathbb{E}_{a_{ji}}\\{a_{ji}\\}\mathbb{E}_{\rho_{m^{\prime}k^{\prime}}}\\{\log\rho_{m^{\prime}k^{\prime}}\\}+\mathbb{E}_{a_{ji}}\\{1-a_{ji}\\}\mathbb{E}_{\rho_{m^{\prime}k^{\prime}}}\\{\log(1-\rho_{m^{\prime}k^{\prime}})\\}\bigg{]}\bigg{(}1-\mathbb{I}_{\\{k^{\prime}\\}}(m^{\prime})\mathbb{I}_{\\{i\\}}(j)\bigg{)}+$
(B.41)
$\displaystyle\qquad\qquad\mathbb{E}_{a_{ii}}\\{a_{ii}\\}\mathbb{E}_{\rho_{k^{\prime}k^{\prime}}}\\{\log\rho_{k^{\prime}k^{\prime}}\\}+\mathbb{E}_{a_{ii}}\\{(1-a_{ii})\\}\mathbb{E}_{\rho_{k^{\prime}k^{\prime}}}\\{\log(1-\rho_{k^{\prime}k^{\prime}})\\}+\mathbb{E}_{\mu_{k^{\prime}}}\left\\{\log\mu_{k^{\prime}}\right\\}\Bigg{\\}}\
+\text{cst.}$ (B.42)

Taking the exponential, and using the distributions derived for
$\hat{q}^{(r-1)}(\rho_{k^{\prime}m^{\prime}})$ and the assumed form of
$\hat{q}^{(r-1)}(a_{ij})$, it follows that the optimal CAVI choice of
$\hat{q}^{(r)}(z_{i}^{\prime})$ is
$\hat{q}^{(r)}(z_{i}^{\prime})=\mathrm{Categorical}(z_{i}^{\prime};\nu_{i})$,
where

$\displaystyle\nu^{(r)}_{ik^{\prime}}\propto\exp\Bigg{\\{}\sum_{j\in\mathcal{V}}\sum_{m^{\prime}\in\mathcal{K}^{\prime}}\nu^{(r-1)}_{jm^{\prime}}\bigg{[}\sigma^{(r-1)}_{ij}\bigg{(}\psi\left(\eta^{(r-1)}_{k^{\prime}m^{\prime}}\right)-\psi\left(\eta^{(r-1)}_{k^{\prime}m^{\prime}}+\zeta^{(r-1)}_{k^{\prime}m^{\prime}}\right)\bigg{)}\
+$ (B.43)
$\displaystyle\quad\left(1-\sigma^{(r-1)}_{ij}\right)\bigg{(}\psi\left(\zeta^{(r-1)}_{k^{\prime}m^{\prime}}\right)-\psi\left(\eta^{(r-1)}_{k^{\prime}m^{\prime}}+\zeta^{(r-1)}_{k^{\prime}m^{\prime}}\right)\bigg{)}+\sigma^{(r-1)}_{ji}\bigg{(}\psi\left(\eta^{(r-1)}_{m^{\prime}k^{\prime}}\right)-\psi\left(\eta^{(r-1)}_{m^{\prime}k^{\prime}}+\zeta^{(r-1)}_{m^{\prime}k^{\prime}}\right)\bigg{)}\
+$ (B.44)
$\displaystyle\quad\left(1-\sigma^{(r-1)}_{ji}\right)\bigg{(}\psi\left(\zeta^{(r-1)}_{m^{\prime}k^{\prime}}\right)-\psi\left(\eta^{(r-1)}_{m^{\prime}k^{\prime}}+\zeta^{(r-1)}_{m^{\prime}k^{\prime}}\right)\bigg{)}\bigg{]}\bigg{(}1-\mathbb{I}_{\\{k^{\prime}\\}}(m^{\prime})\mathbb{I}_{\\{i\\}}(j)\bigg{)}\
+$ (B.45)
$\displaystyle\quad\sigma^{(r-1)}_{ii}\bigg{(}\psi\left(\eta^{(r-1)}_{k^{\prime}k^{\prime}}\right)-\psi\left(\eta^{(r-1)}_{k^{\prime}k^{\prime}}+\zeta^{(r-1)}_{k^{\prime}k^{\prime}}\right)\bigg{)}+\left(1-\sigma^{(r-1)}_{ii}\right)\bigg{(}\psi\left(\zeta^{(r-1)}_{k^{\prime}k^{\prime}}\right)-\psi\left(\eta^{(r-1)}_{k^{\prime}k^{\prime}}+\zeta^{(r-1)}_{k^{\prime}k^{\prime}}\right)\bigg{)}\
+$ (B.46)
$\displaystyle\quad\psi\left(\xi^{(r-1)}_{k^{\prime}}\right)-\psi\left(\sum_{\ell}\xi^{(r-1)}_{\ell}\right)\Bigg{\\}},$
(B.47)

for $i\in\mathbb{V}$ and $k^{\prime}\in\mathcal{K}^{\prime}$.

##  Appendix C Derivation of CAVI updates for the Dirichlet process prior

Under our model specification, the term for $\log p(z_{i}|u)$ can be rewritten
using indicators as

$\displaystyle\log p(z_{i}|u)$
$\displaystyle=\sum_{\ell=1}^{\infty}\bigg{[}\mathbb{I}\\{z_{i}=\ell\\}\log(u_{\ell})+\mathbb{I}\\{z_{i}>\ell\\}\log(1-u_{\ell})\bigg{]}.$
(C.1)

The joint loglikelihood then becomes

$\displaystyle\log p(x,z,u,\lambda)$
$\displaystyle=(\nu-1)\sum_{\ell=1}^{\infty}\log(1-u_{\ell})+(\alpha_{km}-1)\sum_{k,m=1}^{\infty}\bigg{[}\log(\lambda_{km})-\beta_{km}\lambda_{km}\bigg{]}\
+$ (C.2)
$\displaystyle\sum_{i\in\mathcal{V}}\sum_{\ell=1}^{\infty}\bigg{[}\mathbb{I}\\{z_{i}=\ell\\}\log(u_{\ell})\
+\mathbb{I}\\{z_{i}>\ell\\}\log(1-u_{\ell})\bigg{]}\ +$ (C.3)
$\displaystyle\sum_{(i,j)\in\mathcal{E}}\sum_{k,m=1}^{\infty}\tilde{z}_{ik}\tilde{z}_{jm}\bigg{[}x_{ij}(I_{r})\log(\lambda_{km})-\Delta\lambda_{km}\bigg{]}\
+\text{cst.}.$ (C.4)

###  C.1 Approximation of $q(z_{i})$

Recalling that $q(u_{L}=1)=1$, it follows that $\mathbb{I}\\{z_{i}>L\\}=0$. We
can thus compute the necessary expectation as:

$\displaystyle\mathbb{E}_{-z_{i}}$ $\displaystyle\\{\log
p(x,z,u,\lambda)\\}=\mathbb{E}_{-z_{i}}\bigg{\\{}\delta_{z}\sum_{\ell=1}^{L}\bigg{[}\mathbb{I}\\{z_{i}=\ell\\}\log(u_{\ell})+\mathbb{I}\\{z_{i}>\ell\\}\log(1-u_{\ell})\bigg{]}+$
(C.5)
$\displaystyle\quad\sum_{j:(i,j)\in\mathcal{E}}\sum_{k,m=1}^{L}\tilde{z}_{ik}\tilde{z}_{jm}\bigg{(}x_{ij}(I_{r})\log(\lambda_{km})-\Delta\lambda_{km}\bigg{)}\
+$ (C.6)
$\displaystyle\quad\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\sum_{k,m=1}^{L}\tilde{z}_{j^{\prime}k}\tilde{z}_{im}\bigg{(}x_{j^{\prime}i}(I_{r})\log(\lambda_{km})-\Delta\lambda_{km}\bigg{)}-\sum_{k=1}^{L}\tilde{z}_{ik}\bigg{(}x_{ii}(t)\log(\lambda_{kk})-\Delta\lambda_{kk}\bigg{)}\bigg{\\}}+\mathrm{cst.}$
(C.7)
$\displaystyle=\sum_{k=1}^{L}\tilde{z}_{ik}\Bigg{\\{}\delta_{z}\left[\mathbb{E}_{u_{k}}\\{\log(u_{k})\\}+\sum_{\ell=1}^{k-1}\mathbb{E}_{u_{\ell}}\\{\log(1-u_{\ell})\\}\right]+x_{ii}(t)\mathbb{E}_{\lambda_{kk}}\\{\log(\lambda_{kk})\\}-\Delta\mathbb{E}_{\lambda_{kk}}\\{\lambda_{kk}\\}+$
(C.8)
$\displaystyle\quad\sum_{m=1}^{L}\Bigg{[}\sum_{j:(i,j)\in\mathcal{E}}\mathbb{E}_{z_{j}}\\{\tilde{z}_{jm}\\}\bigg{(}x_{ij}(I_{r})\mathbb{E}_{\lambda_{km}}\\{\log(\lambda_{km})\\}-\Delta\mathbb{E}_{\lambda_{km}}\\{\lambda_{km}\\}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}\
+$ (C.9)
$\displaystyle\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\mathbb{E}\\{\tilde{z}_{j^{\prime}m}\\}\bigg{(}x_{j^{\prime}i}(I_{r})\mathbb{E}_{\lambda_{mk}}\\{\log(\lambda_{mk})\\}-\Delta\mathbb{E}_{\lambda_{mk}}\\{\lambda_{mk}\\}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j^{\prime})\bigg{)}\Bigg{]}\Bigg{\\}}+\mathrm{cst.}$
(C.10)

Using the distributions derived for $\hat{q}^{(r-1)}(\lambda_{km})$ and
$\hat{q}^{(r-1)}(u_{\ell})$, it follows that the optimal CAVI distribution is
$\hat{q}^{(r)}(z_{i})=\mathrm{Categorical}\left(z_{i};\tau^{(r)}_{i}\right)$,
with

$\displaystyle\tau^{(r)}_{ik}$
$\displaystyle\propto\exp\Bigg{\\{}\delta_{z}\left[\psi\left(\omega^{(r-1)}_{k}\right)-\psi\left(\omega^{(r-1)}_{k}+\nu^{(r-1)}_{k}\right)+\sum_{\ell=1}^{k-1}\bigg{(}\psi\left(\nu^{(r-1)}_{\ell}\right)-\psi\left(\omega^{(r-1)}_{\ell}+\nu^{(r-1)}_{\ell}\right)\bigg{)}\right]\
+$ (C.11) $\displaystyle\quad
x_{ii}(I_{r})\bigg{(}\psi\left(\alpha^{(r-1)}_{kk}\right)-\log\left(\beta^{(r-1)}_{kk}\right)\bigg{)}-\Delta\frac{\alpha^{(r-1)}_{kk}}{\beta^{(r-1)}_{kk}}+$
(C.12)
$\displaystyle\sum_{m=1}^{L}\Bigg{[}\sum_{j:(i,j)\in\mathcal{E}}\tau^{(r-1)}_{jm}\bigg{(}x_{ij}(I_{r})\left(\psi\left(\alpha^{(r-1)}_{km}\right)-\log\left(\beta^{(r-1)}_{km}\right)\right)-\Delta\frac{\alpha^{(r-1)}_{km}}{\beta^{(r-1)}_{km}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j)\bigg{)}\
+$ (C.13)
$\displaystyle\sum_{j^{\prime}:(j^{\prime},i)\in\mathcal{E}}\tau^{(r-1)}_{j^{\prime}m}\bigg{(}x_{j^{\prime}i}(I_{r})\left(\psi\left(\alpha^{(r-1)}_{mk}\right)-\log\left(\beta^{(r-1)}_{mk}\right)\right)-\Delta\frac{\alpha^{(r-1)}_{mk}}{\beta^{(r-1)}_{mk}}\bigg{)}\bigg{(}1-\mathbb{I}_{\\{k\\}}(m)\mathbb{I}_{\\{i\\}}(j^{\prime})\bigg{)}\Bigg{]}\Bigg{\\}},$
(C.14)

for $i\in\mathcal{V}$ and $k\in\mathcal{K}$.

###  C.2 Approximation of $q(\lambda_{km})$

Recall that $q(u_{L}=1)=1$, and so $\pi(u_{L+1})=\pi(u_{L+2})=\dots$. This
truncation removes the need for an infinite sum. In this way, we can compute
the expectation for the CAVI approximation as:

$\displaystyle\mathbb{E}_{-\lambda_{km}}\\{\log p(x,z,u,\lambda)\\}$
$\displaystyle=\sum_{k,m=1}^{\infty}\mathbb{E}_{-\lambda_{km}}\bigg{[}\delta_{\lambda}\left(\alpha^{(r-1)}_{km}-1\right)\log\lambda_{km}-\delta_{\lambda}\beta^{(r-1)}_{km}\lambda_{km}\
+$ (C.15)
$\displaystyle\quad\sum_{(i,j)\in\mathcal{E}}\tilde{z}_{ik}\tilde{z}_{jm}\bigg{\\{}x_{ij}(I_{r})\log(\lambda_{km})-\Delta\lambda_{km}\bigg{\\}}\bigg{]}+\text{cst.}$
(C.16)
$\displaystyle=\sum_{k,m=1}^{L}\bigg{[}\bigg{(}\delta_{\lambda}\left(\alpha^{(r-1)}_{km}-1\right)+\sum_{(i,j)\in\mathcal{E}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{j}}\\{\tilde{z}_{jm}\\}x_{ij}(I_{r})\bigg{)}\log(\lambda_{km})\
-$ (C.17)
$\displaystyle\quad\bigg{(}\delta_{\lambda}\beta^{(r-1)}_{km}+\Delta\sum_{(i,j)\in\mathcal{E}}\mathbb{E}_{z_{i}}\\{\tilde{z}_{ik}\\}\mathbb{E}_{z_{j}}\\{\tilde{z}_{jm}\\}\bigg{)}\lambda_{km}\bigg{]}+\text{cst.}$
(C.18)

Using the derived distributions for $\hat{q}^{(r-1)}(z_{i})$, it follows that
the optimal CAVI approximation is
$q^{(r-1)}(\lambda_{km})=\mathrm{Gamma}\left(\lambda_{km};\alpha_{km}^{(r)},\beta_{km}^{(r)}\right)$,
where

$\displaystyle\alpha_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\left(\alpha^{(r-1)}_{km}-1\right)+\sum_{(i,j)\in\mathcal{E}}\tau^{(r-1)}_{ik}\tau^{(r-1)}_{jm}x_{ij}(I_{r})+1,$
(C.19) $\displaystyle\beta_{km}^{(r)}$
$\displaystyle=\delta_{\lambda}\beta_{km}^{(r-1)}+\Delta\sum_{(i,j)\in\mathcal{E}}\tau_{ik}^{(r-1)}\tau_{jm}^{(r-1)},$
(C.20)

for $k,m\in\mathcal{K}$.

###  C.3 Approximation of $q(u_{i})$

Again, using that $q^{(r-1)}(u_{L}=1)=1$, we have $\log(1-u_{L})=0$ and
$\mathbb{I}_{\\{k\\}}\\{z_{j}\\}=0$, for $k>L$, the required expectation
becomes:

$\displaystyle\mathbb{E}_{-u_{j}}$ $\displaystyle\\{\log
p(x,z,u,\lambda)\\}=\mathbb{E}_{-u_{j}}\bigg{\\{}\mathbb{I}\\{j<L\\}\bigg{[}\delta_{u}\left(\omega_{j}^{(r-1)}-1\right)\log(u_{j})+\delta_{u}\left(\nu_{j}^{(r-1)}-1\right)\log(1-u_{j})\
+$ (C.21)
$\displaystyle\quad\delta_{z}\sum_{i\in\mathcal{V}}\bigg{(}\mathbb{I}_{\\{j\\}}\\{z_{i}\\}\log(u_{j})+\sum_{k>j}\mathbb{I}_{\\{k\\}}\\{z_{i}\\}\log(1-u_{j})\bigg{)}\bigg{\\}}+\text{cst.}.$
(C.22)
$\displaystyle=\mathbb{I}\\{j<L\\}\bigg{[}\log(u_{j})\bigg{(}\delta_{z}\sum_{i\in\mathcal{V}}\tau^{(r-1)}_{ij}+\delta_{u}\left(\omega_{j}^{(r-1)}-1\right)\bigg{)}\
+$ (C.23)
$\displaystyle\quad\log(1-u_{j})\bigg{(}\delta_{z}\sum_{i\in\mathcal{V}}\sum_{k=j+1}^{L}\tau^{(r-1)}_{ik}+\delta_{u}\left(\nu^{(r-1)}-1\right)\bigg{)}\bigg{]}\
+\mathrm{cst.}$ (C.24)

where we have used the derived forms for $\hat{q}^{(r-1)}(z_{i})$ to compute
the indicator expectations. It follows that the optimal CAVI approximation is
$\hat{q}^{(r)}(u_{j})=\mathrm{Beta}\left(u_{j};\omega_{j}^{(r)},\nu_{j}^{(r)}\right)$,
where

$\displaystyle\omega_{j}^{\prime}$
$\displaystyle=\delta_{z}\sum_{i\in\mathcal{V}}\tau_{ij}^{(r-1)}+\delta_{u}\left(\omega_{j}^{(r-1)}-1\right)+1,$
(C.25) $\displaystyle\nu_{j}^{\prime}$
$\displaystyle=\delta_{z}\sum_{i\in\mathcal{V}}\sum_{k=j+1}^{L}\tau_{ik}^{(r-1)}+\delta_{u}\left(\nu_{j}^{(r-1)}-1\right)+1,$
(C.26)

for $j\in\\{1,\dots,L\\}$.