# Dynamic modeling of spike count data with Conway-Maxwell Poisson variability

Ganchao Wei111Corresponding author<EMAIL_ADDRESS>
Department of Statistics, University of Connecticut  Ian H. Stevenson
Department of Psychological Sciences, University of Connecticut
Department of Biomedical Engineering, University of Connecticut
Connecticut Institute for Brain and Cognitive Science, University of
Connecticut

###### Abstract

In many areas of the brain, neural spiking activity covaries with features of
the external world, such as sensory stimuli or an animal’s movement.
Experimental findings suggest that the variability of neural activity changes
over time and may provide information about the external world beyond the
information provided by the average neural activity. To flexibly track time-
varying neural response properties, here we developed a dynamic model with
Conway-Maxwell Poisson (CMP) observations. The CMP distribution can flexibly
describe firing patterns that are both under- and over-dispersed relative to
the Poisson distribution. Here we track parameters of the CMP distribution as
they vary over time. Using simulations, we show that a normal approximation
can accurately track dynamics in state vectors for both the centering and
shape parameters ($\lambda$ and $\nu$). We then fit our model to neural data
from neurons in primary visual cortex and “place cells” in the hippocampus. We
find that this method out-performs previous dynamic models based on the
Poisson distribution. The dynamic CMP model provides a flexible framework for
tracking time-varying non-Poisson count data and may also have applications
beyond neuroscience.

Keywords— Neural spikes, Non-Poisson count data, Dynamic Conway-Maxwell
Poisson, Laplace approximation

##  1 Introduction

Although many models of neural activity assume that neurons respond with
stable responses to external sensory stimuli or movements, there is
substantial evidence that neural spiking activity changes over time due to
adaptation and plasticity ([1, 2]) as well as spontaneously ([3, 4]). At the
same time, a neuron’s spiking responses on individual trials can be highly
variable, even in the controlled settings with constant stimuli. In most
previous research, trial-to-trial neural variability is assumed to be Poisson
distributed. However, spike count distributions can be substantially more or
less variable than Poisson ([5, 6, 7, 8]), and that the variability also
appears to change over time, in many cases ([9, 10]). Here we introduce a
dynamic model with Conway-Maxwell Poisson observations that can describe non-
Poisson spike statistics and track changing response properties.

Variability appears to be an increasingly important feature of neural
responses and can act as a signature of decision making ([10]), movement
preparation ([11]), or stimulus onset ([9]). Although systems neuroscience has
a long history of studying how external variables influence mean firing rates,
less is known about response variability. Neural activity changes on different
timescales, and distinguishing changes in variability from changes in the mean
response based on sparse, spike observations is a nontrivial statistical
challenge ([12]). Statistical tools to accurately track the sources of
variability within a given experiment may be useful for understanding neural
systems. There has been substantial work developing dynamic Poisson models
([1, 13]), as well as, other Poisson models that can account for fluctuating
response properties with latent or observed variables ([14, 15]). Several
models of neural activity with non-Poisson observations have also been
described ([7, 16, 17]), including a static model with Conway-Maxwell Poisson
observations ([18]). Each of these models, however, whether static or dynamic,
assumes a fixed mean-variance relationship (i.e. fixed dispersion parameters).
Here, to flexibly track how neural variability might change over time, we
explicitly consider changes in both the mean and dispersion.

Here we develop a dynamic GLM with Conway-Maxwell-Poisson (CMP) observations.
The CMP distribution can account for both over- and under-dispersion in spike
count data. To get the closed-form posterior for state vectors with CMP
likelihood, we fit the model using a global Gaussian approximation (Laplace
approximation). Since the state-space of the dynamic model has Markovian
structure, inference is efficient with this approximation, and we estimate the
process noise by maximizing the predictive likelihood. After illustrating the
proposed method in simulations, we apply it to neural activity from primary
visual cortex and place cells in the hippocampus. The dynamic CMP model can
track changes in both the mean and variance of neural responses and
outperforms previous Poisson models.

##  2 Methods

Here we consider a dynamic GLM with Conway-Maxwell Poisson (CMP) observations
to describe time-varying spike counts. We first introduce the model. Although
the CMP distribution allows us to flexibly model non-Poisson variability, one
major challenge with using this model is that there are no closed-form
posteriors for the CMP likelihood. Here, we fit the model using a global
Gaussian approximation, and we discuss several additional technical challenges
that arise when using the CMP distribution with a dynamic GLM. Code is
available at https://github.com/weigcdsb/COM_POISSON.

###  2.1 Dynamic Conway-Maxwell Poisson Model

A count observation $y$, such as the spike count for a neuron, is assumed to
follow the CMP distribution, with parameters $\lambda$ and $\nu$. The
probability mass function (pmf) of CMP distribution is:

$P(Y=y|\lambda,\nu)=\frac{\lambda^{y}}{(y!)^{\nu}}\cdot\frac{1}{Z(\lambda,\nu)}$
(1)

, where $Z(\lambda,\nu)=\sum_{k=0}^{\infty}\frac{\lambda^{k}}{(k!)^{\nu}}$ is
the normalizing constant. The shape parameter $\nu\geq 0$ controls different
dispersion patterns, i.e. equi- $(\nu=1)$, over- $(0\leq\nu<1)$ or under-
dispersion $(\nu>1)$. Three common distributions occur as special cases: 1)
the Poisson $(\nu=1)$, 2) the geometric $(\nu=0,\lambda<1)$, and 3) the
Bernoulli $(\nu\rightarrow\infty)$.

For multiple observations up to T steps, such as simultaneous spike counts
from n neurons, denote the counts at time bin $t$ as
$\bm{y}_{t}=(y_{1t},\ldots,y_{nt})^{\prime}$, for $t=1,\ldots,T$. The
corresponding CMP parameters at $t$ are
$\bm{\lambda}_{t}=(\lambda_{1t},\ldots,\lambda_{nt})^{\prime}$ and
$\bm{\nu}_{t}=(\nu_{1t},\ldots,\nu_{nt})^{\prime}$. Previous work has examined
the CMP-GLM ([19, 20]), and here we focus on the dynamic version of this GLM.
The CMP parameters at $t$ are modeled by two log-linear models,
$\log\bm{\lambda}_{t}=\bm{X}_{t}\bm{\beta}_{t}$ and
$\log\bm{\nu}_{t}=\bm{G}_{t}\bm{\gamma}_{t}$, with
$\bm{\beta}_{t}\in\mathbb{R}^{p}$ and $\bm{\gamma}_{t}\in\mathbb{R}^{q}$, and
$\bm{X}_{t}$ and $\bm{G}_{t}$ denote known predictors. Under the CMP-GLM, the
parameters are static. Here, we assume that they progress linearly with a
Gaussian noise.

The observations follow conditionally independent CMP distributions, given the
state vector
$\bm{\theta}_{t}=(\bm{\beta}^{\prime}_{t},\bm{\gamma}^{\prime}_{t})^{\prime}$.

$\displaystyle\bm{y}_{t}$ $\displaystyle\sim
CMP(\bm{\lambda}_{t},\bm{\nu}_{t})$ (2) $\displaystyle\log\bm{\lambda}_{t}$
$\displaystyle=\bm{X}_{t}\bm{\beta}_{t},\quad\log\bm{\nu}_{t}=\bm{G}_{t}\bm{\gamma}_{t}$

While the state vector $\bm{\theta}_{t}$ evolves linearly with Gaussian noise:

$\displaystyle\bm{\theta}_{1}$ $\displaystyle\sim
N_{p+q}(\bm{\theta}_{0},\bm{Q}_{0})$ (3)
$\displaystyle\bm{\theta}_{t}|\bm{\theta}_{t-1}$ $\displaystyle\sim
N_{p+q}(\bm{F\theta}_{t-1},\bm{Q})$

Given the initial state mean $\bm{\theta}_{0}$, covariance $\bm{Q}_{0}$,
linear dynamics $\bm{F}$ and process covariance $\bm{Q}$.

####  2.1.1 Inference by Gaussian approximation

To fit the model to data we need to estimate the time-varying state vector
$\bm{\Theta}=(\bm{\theta}^{\prime}_{1},\ldots,\bm{\theta}^{\prime}_{T})\in\mathbb{R}^{(p+q)T}$.
In this section, we first assume $\bm{F}$ and $\bm{Q}$ are known. Since the
observations are CMP distributed, we cannot estimate $\bm{\Theta}$ in closed
form. Instead, here we approximate it by a multivariate Gaussian distribution,
$P(\bm{\Theta|Y})\approx N_{(p+q)T}(\bm{\Theta|,\mu,\Sigma})$, with
$\bm{Y}=(\bm{y}^{\prime}_{1},\ldots,\bm{y}^{\prime}_{T})^{\prime}$. The
parameters of this Gaussian are found by a global Laplace approximation, i.e.
$\bm{\mu}=\operatorname*{arg\,max}_{\bm{\Theta}}P(\bm{\Theta|Y})$ and
$\bm{\Sigma}=-(\nabla\nabla_{\bm{\Theta}}\log
P(\bm{\Theta|Y})|_{\bm{\Theta-\mu}})^{-1}$. The log-posterior is given by:

$\displaystyle\log P(\bm{\Theta|Y})$
$\displaystyle=\sum_{t=1}^{T}l_{t}-\frac{1}{2}(\bm{\theta}_{1}-\bm{\theta}_{0})^{\prime}\bm{Q}_{0}^{-1}(\bm{\theta}_{1}-\bm{\theta}_{0})-\frac{1}{2}\sum_{t=2}^{T}(\bm{\theta}_{t}-\bm{F\theta}_{t-1})^{\prime}\bm{Q}^{-1}(\bm{\theta}_{t}-\bm{F\theta}_{t-1})$
(4) $\displaystyle l_{t}$ $\displaystyle=l(\bm{\theta}_{t})=\log
P(\bm{y}_{t}|\bm{\theta}_{t})=\sum_{i=1}^{n}y_{it}\log\lambda_{it}-\nu_{it}y_{it}!-\log
Z(\lambda_{it},\nu_{it})$

, where $l(\cdot)$is the log-likelihood. The log-posterior is concave ([21]),
and the Markovian structure of the state vector dynamics makes it possible to
optimize by Newton-Raphson (NR) in $\mathcal{O}(T)$ time ([22]).After the
Newton update, we can further quantify the uncertainty for the CMP parameters
and the underlying rates, as in Appendix A.

There are several technical challenges involved with performing the Newton
update with CMP observations. Firstly, in order to find the gradient and
Hessian we need to calculate moments of $Y_{it}$ and $\log Y_{it}!$, which
have no closed forms ([23]). We can calculate these moments by truncated
summation. However, when $\lambda\geq 2$ and $\nu\leq 1$, truncated summation
is computationally costly since we need many steps for accurate approximation.
In this case, we approximate the moments using previous ([19, 24]) asymptotic
results as in Appendix B. A second challenge is that the Hessian is not robust
to outliers. Outliers often result in the Hessian being close to singular or
even positive-definite. See details in Appendix C. To ensure robustness, we
use Fisher scoring where the observed information is replaced by the expected
information. Finally, a third challenge is that the Newton updates take a long
time to converge if the initial state estimate is far from the maximum of the
posterior, especially when $T$ is large. To resolve this issue, we use a
smoothing estimate with local Gaussian approximation as a “warm start”.
Forward filtering for a dynamic Poisson model has been previously described in
[13], and here we implement CMP filtering following the same rationale. Let
$\bm{\theta}_{t|t-1}=E(\bm{\theta}_{t}|\bm{y}_{1},\ldots,\bm{y}_{t-1})$ and
$\bm{\Sigma}_{t|t-1}=Var(\bm{\theta}_{t}|\bm{y}_{1},\ldots,\bm{y}_{t-1})$ be
the mean and variance for the one-step prediction density and
$\bm{\theta}_{t|t}=E(\bm{\theta}_{t}|\bm{y}_{1},\ldots,\bm{y}_{t})$ and
$\bm{\Sigma}_{t|t}=Var(\bm{\theta}_{t}|\bm{y}_{1},\ldots,\bm{y}_{t})$ be mean
and variance for the posterior density, then the filtering update for step $t$
is given by

$\displaystyle\bm{\theta}_{t|t-1}$ $\displaystyle=\bm{F\theta}_{t-1|t-1}$ (5)
$\displaystyle\bm{\Sigma}_{t|t-1}$
$\displaystyle=\bm{F\Sigma}_{t-1|t-1}\bm{F}^{\prime}+\bm{Q}$
$\displaystyle\bm{\theta}_{t|t}$
$\displaystyle=\bm{\theta}_{t|t-1}+(\bm{\Sigma}_{t|t})\left[\frac{\partial
l_{t}}{\partial\bm{\theta}_{t}}\right]_{\bm{\theta}_{t|t-1}}$
$\displaystyle(\bm{\Sigma}_{t|t})^{-1}$
$\displaystyle=(\bm{\Sigma}_{t|t-1})^{-1}-\left[\frac{\partial^{2}l_{t}}{\partial\bm{\theta}_{t}\partial\bm{\theta}^{\prime}_{t}}\right]_{\bm{\theta}_{t|t-1}}$

Here, to again ensure robustness, we use Fisher scoring when updating the
state covariance. We then find smoothed estimates using a backward pass
([25]). Although doing smoothing is fast, the estimates can be inaccurate,
especially when there are large changes in the state vector. In the forward
filtering stage, the Gaussian approximation at each step $t$ is conducted
locally at the recursive prior $\bm{\theta}_{t|t-1}$. This will be
statistically inefficient when the recursive prior is too far away from the
posterior mode, or when there is a large change in the state vector. Moreover,
Fisher scoring reduces the efficiency of the smoother even further. The
smoother provides reasonable initial estimates, but estimation accuracy is
substantially improved by using Newton’s method to find the global Laplace
approximation for the posterior.

####  2.1.2 Estimating process noise

For the applications to neural data examined here, we assume that
$\bm{F}=\bm{I}$. However, we still need to estimate the process noise
$\bm{Q}$. When $n$ is small, especially when $n=1$, different $\bm{Q}$ values
will have a substantial influence on estimation. One possible way to estimate.
One possible way to estimate $\bm{Q}$ is to use an Expectation Maximization
(EM) algorithm as in [26]. However, using the Laplace approximation for
$\bm{\Theta}$ during E-step breaks the usual guarantee of non-decreasing
likelihoods in EM, and, hence, may lead to divergence. To avoid that, we could
sample the posterior directly by MCMC. However, the lack of closed-form
moments for the CMP distribution makes sampling computationally intensive.
Here, to estimate $\bm{Q}$ robustly and quickly, we instead assume $\bm{Q}$ is
diagonal and estimate it by maximizing the prediction likelihood in the
filtering stage, as in [27].

##  3 Results

###  3.1 Tracking the mean and dispersion of spike counts over time

To illustrate how the dynamic CMP model can track both time-varying mean and
dispersion, we simulated a neuron with a time-varying tuning curve, where the
response to 100 evenly-spaced hypothetical visual stimuli shifts over 100
trials. Here, the neuron’s tuning curve is determined by a linear combination
of cubic B-spline basis functions with equally-spaced knots. The stimulus that
evokes the highest average response – the “preferred orientation” – is
initially 80 deg, but shifts over the course of the experiment, and the
response amplitude also increases over time (Fig. 1A). Meanwhile, the
dispersion pattern also changes: the responses are initially over-dispersed
relative to a Poisson distribution and then become under-dispersed (Fig. 1B).
Noisy observations are sampled from the Conway-Maxwell Poisson distribution at
each time (Fig. 1C), mimicking the types of experimental observations
collected during adaptation experiments in primary visual cortex ([28]). We
then fitted the simulated spike observations using the same predictor
variables as the generative model: the covariates for $\bm{\lambda_{t}}$
capture the tuning curve with $\bm{X}_{t}$ as a 10-knot cubic spline basis
expansion of the orientation, and the covariate for the shape parameter
$\bm{\nu_{t}}$ does not depend on the stimulus orientation $\bm{G}_{t}$. The
fitted results match the ground truth well, for both the mean (Fig. 1C) and
Fano factors (variance-to-mean ratio, Fig. 1D).

This model-based approach provides estimates of tuning curves and dispersion
at each time point. In cases where the tuning curve and variability change
simultaneously, this approach can efficiently track both. By using the model
with CMP observations, rather than Poisson or negative-binomial observations,
the Fano factor can be both ¡1 (under-dispersed) and ¿1 (over-dispersed).

Figure 1: A simulated neuron with a shifting firing and dispersion pattern.
(A) The tuning curve of the neuron shifts over time, with the preferred
stimulus orientation changing and the response amplitude increasing. (B) At
the same time, the variability in spiking changes from being over-dispersed
relative to a Poisson distribution to under-dispersed. This leads to a
decreasing Fano factor (variance-to-mean ratio) from 1.9 to 0.4 overall. (C)
To illustrate the shifts, we show the tuning curve at two time points: Trial
20 (blue) and Trial 80 (red). Dots denote observed spike counts. The solid
lines are the ground truth in mean firing rate, while the corresponding dashed
lines are the fitted values. (D) To illustrate the shift in dispersion over
time we show the true (solid) and estimated (dashed) Fano factor for two
specific stimuli as a function of time. The dispersion for the early preferred
orientation is shown in cyan, while the dispersion for the late preferred
orientation is shown in yellow.

Changes in tuning have been widely documented in systems neuroscience both due
to changing environment and spontaneous nonstationarity. Changes in
variability also occur, but have been less well studied. With the CMP model,
the mean and dispersion are both tracked and, thus, changes in variability can
occur even when the mean is stable. To illustrate this potential, we simulated
a neuron whose mean firing rate is controlled to be constant, but whose Fano
factor varies over time. Here $\bm{X}_{t}$ is a 5-knot cubic B-spline basis
expansion of the orientation and $\bm{G}_{t}=1$. The model recovers the true
mean firing rate (Fig. 2A) and capture the fluctuations in variance (Fig. 2B)
at the same time. However, the estimated Fano factor is somewhat oversmoothed
when the process noise $\bm{Q}$ is optimized by maximizing the predictive
likelihood (see 2.1.2).

Figure 2: Constant mean with fluctuations in dispersion. (A) The first two
panels show the true mean firing rate and the simulated observations. The last
panel show the fitted mean response. (B) Although the mean response is
constant, the Fano factor varies across the trial (black line). The colored
line show the fitted result.

Although dynamic Poisson models have been applied in some neuroscientific
settings, when spike counts are not Poisson distributed the model estimates
can be biased. Since the dispersion influences estimates of the process noise
$\bm{Q}$, estimates of the mean in the dynamic Poisson model can be effected
by over- or under-dispersion. To illustrate this interaction here we simulate
a place cell from the hippocampus whose “place field” drifts over time. The
true mean is determined by a Gaussian function where the preferred position
varies over time. The spike counts are then generated by CMP distributions,
here over-dispersed with constant shape parameter $\bm{\nu}_{t}=0.1$. We fit
1000 observations randomly sampled from 100 “runs” of a linear track. We find
that, in this data-limited regime, the dynamic Poisson model and the dynamic
CMP model give substantially different estimates of the time-varying place
field (Fig. 3A). The dynamic Poisson model, in this case, under-estimates the
firing rate at the true preferred position and under-estimates the uncertainty
(Fig. 3B).

Figure 3: Dynamic CMP and dynamic Poisson estimates differ. Simulated over-
dispersed place cell spiking is generated by the dynamic CMP model with
$\nu_{t}=0.1$. (A) We then fit dynamic CMP and dynamic Poisson models with
$\bm{X}_{t}$, a 2-knot B-spline expansion for position ($\bm{G}_{t}=1$ for the
CMP). (B) When evaluating the response at the true preferred position for each
run, the dynamic Poisson estimates are biased (under-estimated) and the
uncertainty is also underestimated. The solid line gives the MAP estimates of
mean firing rate, and the dashed lines show one S.D. credible intervals. The
standard deviations of dynamic CMP estimates are calculated using the
truncated summations (see details in APPENDIX B), while the standard
deviations for the dynamic Poisson model are from a log-normal distribution.

###  3.2 Application to Experimental Data

We next applied our method to three publicly available datasets of
extracellular spike recordings: 1) Utah array recordings of visually evoked
activity from anesthetized macaque primary visual cortex (“V1 data”), 2)
multi-shank silicon probe recordings from hippocampus of a rat running back-
and-forth on a linear maze (“HC data”) and 3) a speed-tuned neuron recorded
from the anterior pretectal nucleus (APN) of an awake mouse (”APN neuron”)

####  3.2.1 V1 Data

In the V1 dataset. CRCNS pvc-11 ([29]), anesthetized macaque monkeys viewed
full-field sinusoidal grating movies while neural activity was recorded by a
96-channel “Utah” array. Extracellular spiking activity was recorded on each
electrode, and spike waveform segments were sorted by hand with modified
competitive mixture decomposition methods ([30]). Here we use data from one
animal (Monkey 1) viewing a movie of drifting sinusoidal gratings with 100
different drift directions presented in pseudorandom order (300ms each, 30s
movie in total), and the movie was repeated 120 times. Here we analyze spike
counts following each stimulus presentation from the period 50-350ms after
stimulus onset, considering the response delay. For further details on how the
data were obtained, see [31, 32].

As with many neurons in visual cortex, the responses of the neurons in this
dataset are tuned to the stimulus direction. Neurons respond to some
directions of stimuli more than others, but the spike counts from trial to
trial are highly variable. Here, we are specifically interested in tracking
changes tuning curves and changes in variability over time. Fig. 4A shows
responses from one example neuron with a preferred direction around 240 deg.
This neuron is somewhat direction insensitive, and also responds with
increased spiking to stimuli moving in the opposite direction, around 70 deg.
After fitting the dynamic CMP model to these data, we find that the tuning
curve itself appears mostly stable, but the overall firing rate increases over
the course of the recording (Fig. 4B). At the same time, the Fano factor
decreases over the course of the recording (Fig. 4C).

Although the data here is structured in 120 “trials” the data are collected
sequentially, and we model nonstationary at the level of individual
observations. For the predictors $\bm{X}_{t}$ and $\bm{G}_{t}$ we use cubic
B-spline basis functions with periodic boundary conditions over the grating
directions. Results for the example neuron use 5 and 3 equally-spaced knots
for $\bm{X}_{t}$ and $\bm{G}_{t}$, respectively. Fitting the model with half
of the data (in a speckled hold-out pattern) gives patterns for the mean
response (Fig. 4B) and Fano factor (Fig. 4C) that are similar to those using
the full data. However, since the model-based approach provides a continuous
estimate of the state vectors, the Fano factor estimated by the dynamic CMP
model differs from a simple estimate of the Fano factor calculated using a
sliding window (Fig. 4C).

We then compare the performance of multiple models on data from all 74 neurons
in this recording (Fig. 4D). We assess four dynamic models: (1) dynamic CMP,
with 5 knots for $\bm{X}_{t}$ and 3 knots for $\bm{G}_{t}$, denoted as
dCMP-(5,3); (2) dynamic CMP with $\bm{G}_{t}=1$, dCMP-(5,1); (3) dynamic CMP
with constant $\nu_{t}$, dCMP-(5)-$\nu$ (fit by coordinate descent) and (4) a
dynamic Poisson model, dPoi-(5). Additionally, we assess three static models:
(1) static CMP, sCMP-(5,3); (2) static CMP with $\bm{G}_{t}=1$, sCMP-(5,1) and
(3) static Poisson, sPoi-(5). The held-out log-likelihoods relative to a
homogeneous static Poisson model show that the CMP-based models, both dynamic
and static, outperform the Poisson-based models (Fig. 4D). The dynamic models
perform slightly better than the corresponding static models, on average. The
best performance on test data comes from modeling nonstationarity and
stimulus-dependence with the full dynamic CMP model dCMP-(5,3). However, the
benefit of adding nonstationary shape parameter (dCMP-(5)-$\nu$ vs.
dCMP-(5,1)) and of adding stimulus-dependent shape parameter (dCMP-(5,1) vs
dCMP-(5,3)) tend to be small for these data.

Figure 4: Modeling nonstationary spiking from visual cortex. (A) The spike
counts of one example neuron from V1 in response to drifting grating stimuli
with different drift directions presented over 120 trials. Two preferred
directions (estimated by the CMP model) are marked by the dashed lines. (B)
Estimated mean for the dynamic CMP model dCMP-(5,3) when fit to all the data
and only half of the observations (held out in a speckled pattern). (C) Fano
factor estimates for the two models, along with a direct estimate from
15-trial sliding windows, at the two preferred directions. Dashed lines denote
±1 standard deviation around the window estimates, obtained by Bayesian
bootstrapping. (D) Model comparison for all 74 neurons in the V1 dataset. In
these models, 4 are dynamic and the remaining 3 are static, with different
noise distributions (Poisson vs CMP) and bases. The training and test log-
likelihood ratios (bits/spike) with respect to a homogeneous static Poisson
model are shown for all neurons in grey lines. The solid orange lines and
numbers denote the medians, and the dashed lines show the first and third
quartiles.

####  3.2.2 HC Data

In the HC dataset, CRCNS hc-3 ([33]), a rat was running back and forth along a
250cm linear track. Extra cellular spiking activity was recorded in dorsal
hippocampus using multi-shank silicon probes. Spikes were automatically
sorting using KlustaKwik followed by manual adjustment ([34]). Here we use
data from one 66 min recording session (ec014-468) and analyze spike counts in
200ms bins. For further details on how the data were obtained, see [35]

As with many neurons in hippocampus, the responses of the neurons in this
dataset are tuned to the rat’s position along the track. Neurons spike at
specific locations, but the place fields can also shift over time and the
spike counts from run to run are highly variable. Fig. 5A shows an example
from one neuron with two place fields where the location and firing within the
place field vary over the course of the recording. Compared to the data from
V1, neural responses of place cells in hippocampus tend to be sparser and more
selective. Many place cells also tend to be direction tuned – spiking only
when the animal is running in one direction down the track but not the other.
We, thus, fit the data using a dynamic CMP model with 12 equally-spaced knots
for $\bm{X}_{t}$ with a circular representation of position, and let
$\bm{G}_{t}=1$.

For this example neuron, the dynamic CMP model accurately tracks the time-
varying place field (Fig. 5B). We then evaluate the fitted Fano factors at the
peaks of the two place fields (Fig. 5C). Compared to example from the V1 data,
the spiking of this example place cell is much more highly dispersed. The Fano
factors vary over time and are also specific in both position and running
direction. We then compare model performance on 78 neurons from this recording
(19 neurons were excluded due to sparse spiking patterns). In these data, the
dynamic models are generally better than the static models (Fig. 5D). Within
groups of dynamic or static models, CMP-based models are consistently better
than the Poisson-based models.

Figure 5: Modeling nonstationarity in hippocampal spiking activity. (A) Spike
counts of an example neuron from the hippocampus recorded while a rat was
running back-and-forth on a linear maze. The blue lines show the animal’s
position over time. The black circles denote spike counts with the radius of
each circle corresponding to number of spikes. (B) The heatmap shows fitted
mean firing rate for the dynamic CMP model dCMP-(12,1). The colored lines show
peaks for two place fields, chosen based on the model fit. (C) The estimated
Fano factors at the two place field peaks, with each running direction (+ vs.
-) shown separately. (D) The training and test log-likelihood ratios
(bits/spike) with respect to a homogeneous static Poisson model for 78 neurons
(gray lines). Here the test set log-likelihood ratios are calculated using 5%
of the data held-out in a speckled pattern. The red solid lines and numbers
denote median values, while the dashed lines show the first and third
quartiles.

####  3.2.3 APN Neuron

To illustrate a case where the dynamic CMP provides a qualitatively better
description of neural activity compared to previous models we show results
from one neuron recorded from the Allen Institute Visual Coding Neuropixels
dataset. See detailed data description in [36]. Here, when examining tuning to
running speed we found a neuron in APN, whose responses were speed tuned -
increasing firing with increasing running speed (Fig. 6A), but also highly
under-dispersed relative to a Poisson distribution (Fano factor $<$1). We
analyzed spike counts in 200ms bins for the 160 min recording
(ecephys_session_id 719161530, unit_id 950917034), and fit the spiking
activity for the whole session using the dynamic CMP model. Here we use a
nonlinear function of the running speed $v_{t}$ for the
$\bm{X}_{t}=[1\,f(v_{t})]$ ($f(v)=v/(1+0.1v)$) and $\bm{G}_{t}=[1\,v_{t}]$.

Evaluating all of the data and averaging over time, the neuron is
significantly under-dispersed (Fano factors less than 1), and both the mean
response and dispersion appear tuned to running speed (Fig. 6B). However, we
found that this neuron’s speed tuning is somewhat nonstationary, with the
baseline firing rate shifting over time (Fig. 6C). Within individual segments
of the recording the Fano factor is much lower than the Fano factor evaluated
for the entire recording. Moreover, the Fano factor within each chunk doesn’t
show a strong relationship to running speed (Fig. 6D). This suggests that the
apparent relationship between running speed and the overall Fano factor is a
byproduct of the nonstationarity. Using the static CMP model may be able to
describe the underdispersion of these responses, but would miss this key
feature of the data by assuming that the tuning curve is static. The dynamic
Poisson model, on the other hand, is able to describe the nonstationary in the
tuning curve, but cannot describe the underdispersion of these data, since the
Poisson model always assumes Fano factor = 1.

Figure 6: Modeling a speed-tuned neuron in mouse anterior pretectal nucleus
with the dynamic CMP model Here we analyze a 160 min recording of a neuron
from APN, with spike counts calculated in 200ms bins. The whole session was
evenly partitioned into chunks for visualization. (A) Spike counts and mouse
running speed for the recording session. (B) The spiking activity and Fano
factor as a function of speed for the whole session. Dots denote results from
observed spike counts. Errorbars denote 95% confidence intervals from Bayesian
bootstrapping. Red lines denote fits from the dynamic CMP model averaged over
time. (C) The top row shows 2D histograms for spiking counts vs. running speed
for each chunk. The second row shows corresponding model fitted densities
(averages within each chunk). The lines show predicted average spike counts
(tuning curves) for each chunk. (D) There is substantial variation in the
tuning curves across chunks, and the Fano factor tends to be smaller on
individual chunks than in the overall recording.

##  4 Discussion

Here we introduced a dynamic CMP model to track changes in both the mean and
dispersion of neural spikes over time. A global Laplace approximation with a
smoothing-based initialization can provide accurate and computationally
efficient model estimates. In both simulations and applications with
experimental data we find that this model out-performs previous static and
dynamic Poisson models, and may, thus, be a useful tool for understanding the
role of variability in neural systems. While many studies have characterized
noise ([7, 37, 38]) or non-stationarity ([4, 39]) separately, modeling changes
in the mean-variance relationship directly may allow us to more closely
examine the role of variability in the brain.

The extent to which the dynamic CMP model can predict neural responses more
accurately than the dynamic Poisson model or static non-Poisson models depends
on the neural activity itself. Here with the V1 data we found a 6% improvement
in test log-likelihood ratios between the dynamic and static CMP, while in the
HC data there is a 35% improvement. The extent of spontaneous changes in
neural responses is somewhat unclear, some evidence suggests that neurons can
be relatively stable in some circumstances ([40, 41, 42]), but variability
appears to differ across brain areas ([43]). More accurate spike sorting may
account for some degree of instability ([44]), and the degree of spontaneous
changes may also depends on the brain area ([45]). However, neurons do clearly
change both their average responses and dispersion in many situations.

Although the current model works well for fitting neural spikes, there are
some potential improvements. First, the state transition matrix $\bm{F}$ is
currently assumed to be known and is fixed to $\bm{I}$ for convenience in our
simulations and experimental analysis. This doesn’t allow for interactions
between state vectors and may limit the usage in some situations. However,
when using a Laplace approximation for the state vectors, $\bm{F}$ can be
estimated using the EM algorithm as in [26]. Secondly, although the CMP
distribution can flexibly model over- and under-dispersed data, the assumed
mean-variance relationship may not be appropriate in some cases. To more
flexibly model the dispersion for single observations, it may be useful to
instead consider the generalized count (GC) distribution ([46]), which
includes the CMP distribution as a special case. This model has been applied
in the context of linear dynamical systems, similar to the dynamical factor
analysis model, with a fixed dispersion function ([16]). To track fluctuations
in dispersion more flexibly, it could be useful to allow the function to vary
dynamically similar to $\nu_{t}$ in the dynamic CMP model here.

The best modeling strategy also likely depends on what researchers want to
know about the variability. Omitted variables ([47, 48]) and history effects
([49]) can alter apparent observation noise. For example, in the hippocampus,
place cell firing is highly variable on different passes through the field
([50]). This may be partially due to joint selectivity to position, speed, and
head direction, as well as the influence of local field potentials. Here,
rather than model these distinct covariates assuming Poisson observations, we
allow the variability to be non-Poisson and introduce a dynamic GLM with CMP
observations. However, doubly stochastic Poisson models ([51]) or latent
variable models with fixed mean-variance relationships ([16, 17]) may also be
able to account for some differences in the variance over time. Nonetheless,
the dynamic CMP model may provide a useful tool for neuroscientists to study
the role of variance directly. Since the static CMP model can improve decoding
of external variables in some cases ([52]), the dynamic CMP may lead to
further improvements in decoding by tracking nonstationarity in neural
response properties.

## Acknowledgments

This material is based upon work supported by the National Science Foundation
under Grant No. 1931249.

## References

  * [1] Emery N. Brown, David P. Nguyen, Loren M. Frank, Matthew A. Wilson, and Victor Solo.  An analysis of neural receptive field plasticity by point process adaptive filtering.  Proceedings of the National Academy of Sciences, 98(21):12261–12266, oct 2001.
  * [2] Nicholas A. Lesica, Jianzhong Jin, Chong Weng, Chun I. Yeh, Daniel A. Butts, Garrett B. Stanley, and Jose Manuel Alonso.  Adaptation to Stimulus Contrast and Correlations during Natural Visual Stimulation.  Neuron, 55(3):479–491, 2007.
  * [3] U Rokni, A G Richardson, E Bizzi, and H S Seung.  Motor learning with unstable neural representations.  Neuron, 54(4):653–666, 2007.
  * [4] George J Tomko and Donald R Crapper.  Neuronal variability: non-stationary responses to identical visual stimuli.  Brain Research, 79(3):405–418, 1974.
  * [5] Gaby Maimon and John a. Assad.  Beyond Poisson: Increased Spike-Time Regularity across Primate Parietal Cortex.  Neuron, 62(3):426–440, 2009.
  * [6] Asohan Amarasingham, Ting Li Chen, Stuart Geman, Matthew T. Harrison, and David L. Sheinberg.  Spike Count Reliability and the Poisson Hypothesis.  Journal of Neuroscience, 26(3):801–809, jan 2006.
  * [7] Michael R DeWeese, Michael Wehr, and Anthony M Zador.  Binary spiking in auditory cortex.  The Journal of neuroscience : the official journal of the Society for Neuroscience, 23(21):7940–7949, 2003.
  * [8] P Kara, P Reinagel, and R C Reid.  Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons.  Neuron, 27(3):635–646, 2000.
  * [9] Mark M Churchland, Byron M Yu, John P Cunningham, Leo P Sugrue, Marlene R Cohen, Greg S Corrado, William T Newsome, Andrew M Clark, Paymon Hosseini, Benjamin B Scott, David C Bradley, Matthew a Smith, Adam Kohn, J Anthony Movshon, Katherine M Armstrong, Tirin Moore, Steve W Chang, Lawrence H Snyder, Stephen G Lisberger, Nicholas J Priebe, Ian M Finn, David Ferster, Stephen I Ryu, Gopal Santhanam, Maneesh Sahani, and Krishna V Shenoy.  Stimulus onset quenches neural variability: A widespread cortical phenomenon.  Nature Neuroscience, 13(3):369–378, 2010.
  * [10] Anne K Churchland, R. Kiani, R. Chaudhuri, Xiao Jing Wang, Alexandre Pouget, and M. N. Shadlen.  Variance as a Signature of Neural Computations during Decision Making.  Neuron, 69(4):818–831, 2011.
  * [11] M M Churchland, B M Yu, S I Ryu, G Santhanam, and K V Shenoy.  Neural Variability in Premotor Cortex Provides a Signature of Motor Preparation.  Journal of Neuroscience, 26(14):3697, 2006.
  * [12] Michael DeWeese and Anthony Zador.  Asymmetric Dynamics in Optimal Variance Adaptation.  Neural Computation, 10(5):1179–1202, 1998.
  * [13] Uri T Eden, Loren M. Frank, Riccardo Barbieri, Victor Solo, and Emery N. Brown.  Dynamic Analysis of Neural Encoding by Point Process Adaptive Filtering.  Neural Computation, 16(5):971–998, may 2004.
  * [14] Gabriela Czanner, Uri T Eden, Sylvia Wirth, Marianna Yanike, Wendy A Suzuki, and Emery N Brown.  Analysis of between-trial and within-trial neural spiking dynamics.  Journal of Neurophysiology, 99(5):2672–2693, 2008.
  * [15] Anne C. Smith and Emery N. Brown.  Estimating a state-space model from point process observations.  Neural Computation, 15(5):965–991, may 2003.
  * [16] Yuanjun Gao, Lars Buesing, Krishna V Shenoy, and John P Cunningham.  High-dimensional neural spike train analysis with generalized count linear dynamical systems.  Technical report, 2015.
  * [17] Jonathan W. Pillow and James G Scott.  Fully Bayesian inference for neural models with negative-binomial spiking.  Technical report, 2012.
  * [18] Ian H. Stevenson.  Flexible models for spike count data with both over- and under- dispersion.  Journal of Computational Neuroscience, 41(1):29–43, aug 2016.
  * [19] Suneel Babu Chatla and Galit Shmueli.  Efficient estimation of COM–Poisson regression and a generalized additive model.  Computational Statistics and Data Analysis, 121:71–88, may 2018\.
  * [20] Kimberly F. Sellers and Galit Shmueli.  A flexible regression model for count data.  The Annals of Applied Statistics, 4(2):943 – 961, 2010.
  * [21] Ramesh C. Gupta, S. Z. Sim, and S. H. Ong.  Analysis of discrete data by Conway–Maxwell Poisson distribution.  AStA Advances in Statistical Analysis, 98(4):327–343, oct 2014\.
  * [22] Liam Paninski, Yashar Ahmadian, Daniel Gil Ferreira, Shinsuke Koyama, Kamiar Rahnama Rad, Michael Vidne, Joshua Vogelstein, and Wei Wu.  A new look at state-space models for neural data.  J. Comput. Neurosci., 29(1–2):107–126, aug 2010.
  * [23] Galit Shmueli, Thomas P. Minka, Joseph B. Kadane, Sharad Borle, and Peter Boatwright.  A useful distribution for fitting discrete data: revival of the Conway–Maxwell–Poisson distribution.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(1):127–142, jan 2005.
  * [24] Robert E Gaunt, Satish Iyengar, Adri B Olde Daalhuis, Burcin Simsek, and B E Robert Gaunt.  An asymptotic expansion for the normalizing constant of the Conway-Maxwell-Poisson distribution.  Ann Inst Stat Math, 71:163–180, 2019.
  * [25] H E Rauch, F Tung, and C T Striebel.  Maximum likelihood estimates of linear dynamic systems.  AIAA Journal, 3(8):1445–1450, aug 1965.
  * [26] Jakob H. Macke, Lars Buesing, John P. Cunningham, Byron M. Yu, Krishna V. Shenoy, and Maneesh Sahani.  Empirical models of spiking in neural populations.  Advances in Neural Information Processing Systems, 24, 2011.
  * [27] Ganchao Wei and Ian H. Stevenson.  Tracking Fast and Slow Changes in Synaptic Weights From Simultaneously Observed Pre- and Postsynaptic Spiking.  Neural computation, 33(10):2682–2709, sep 2021.
  * [28] Valentin Dragoi, Jitendra Sharma, and Mriganka Sur.  Adaptation-Induced Plasticity of Orientation Tuning in Adult Visual Cortex.  Neuron, 28(1):287–298, oct 2000.
  * [29] A. Kohn and M.A. Smith.  Utah array extracellular recordings of spontaneous and visually evoked activity from anesthetized macaque primary visual cortex (V1), 2016.
  * [30] Shy Shoham, Matthew R. Fellows, and Richard A. Normann.  Robust, automatic spike sorting using mixtures of multivariate t-distributions.  Journal of neuroscience methods, 127(2):111–122, aug 2003.
  * [31] Ryan C. Kelly, Matthew A. Smith, Robert E. Kass, and Tai Sing Lee.  Local field potentials indicate network state and account for neuronal response variability.  Journal of Computational Neuroscience, 29(3):567–579, dec 2010\.
  * [32] Matthew A. Smith and Adam Kohn.  Spatial and temporal scales of neuronal correlation in primary visual cortex.  Journal of Neuroscience, 28(48):12591–12603, nov 2008.
  * [33] K. Mizuseki, A. Sirota, E. Pastalkova, K. Diba, and G. Buzsáki.  Multiple single unit recordings from different rat hippocampal and entorhinal regions while the animals were performing multiple behavioral tasks.  CRCNS.org, 2013.
  * [34] Cyrille Rossant, Shabnam N. Kadir, Dan F.M. Goodman, John Schulman, Maximilian L.D. Hunter, Aman B. Saleem, Andres Grosmark, Mariano Belluscio, George H. Denfield, Alexander S. Ecker, Andreas S. Tolias, Samuel Solomon, György Buzski, Matteo Carandini, and Kenneth D. Harris.  Spike sorting for large, dense electrode arrays.  Nature Neuroscience 2016 19:4, 19(4):634–641, mar 2016.
  * [35] Kenji Mizuseki, Kamran Diba, Eva Pastalkova, Jeff Teeters, Anton Sirota, and György Buzsáki.  Neurosharing: large-scale data sets (spike, LFP) recorded from the hippocampal-entorhinal system in behaving rats.  F1000Research 2014 3:98, 3:98, jul 2014.
  * [36] Joshua H. Siegle, Xiaoxuan Jia, Séverine Durand, et al.  Survey of spiking in the mouse visual system reveals functional hierarchy.  Nature 2021 592:7852, 592(7852):86–92, jan 2021.
  * [37] Michael R Deweese and Anthony M Zador.  Shared and private variability in the auditory cortex.  Journal of neurophysiology, 92(3):1840–1855, 2004.
  * [38] Wahiba Taouali, Giacomo Benvenuti, Pascal Wallisch, Frédéric Chavane, and Laurent U Perrinet.  Testing the odds of inherent vs. observed overdispersion in neural spike counts.  Journal of neurophysiology, 115(1):434–44, jan 2016.
  * [39] W Wu and N G Hatsopoulos.  Real-time decoding of nonstationary neural activity in motor cortex.  Neural Systems and Rehabilitation Engineering, IEEE Transactions on, 16(3):213–222, 2008.
  * [40] Cynthia A Chestek, Aaron P Batista, Gopal Santhanam, Byron M Yu, Afsheen Afshar, John P Cunningham, Vikash Gilja, Stephen I Ryu, Mark M Churchland, and Krishna V Shenoy.  Single-neuron stability during repeated reaching in macaque premotor cortex.  Journal of Neuroscience, 27(40):10742–10750, 2007.
  * [41] Ian H. Stevenson, Anil Cherian, Brian M. London, Nicholas A. Sachs, Eric Lindberg, Jacob Reimer, Marc W. Slutzky, Nicholas G. Hatsopoulos, Lee E. Miller, and Konrad P. Kording.  Statistical assessment of the stability of neural movement representations.  Journal of Neurophysiology, 106(2):764–774, aug 2011.
  * [42] Adam S Dickey, Aaron Suminski, Yali Amit, and Nicholas G Hatsopoulos.  Single-unit stability using chronically implanted multielectrode arrays.  Journal of Neurophysiology, 102(2):1331–1339, 2009.
  * [43] Yasuhiro Mochizuki, Tomokatsu Onaga, Hideaki Shimazaki, Takeaki Shimokawa, Yasuhiro Tsubo, Rie Kimura, Akiko Saiki, Yutaka Sakai, Yoshikazu Isomura, Shigeyoshi Fujisawa, Ken Ichi Shibata, Daichi Hirai, Takahiro Furuta, Takeshi Kaneko, Susumu Takahashi, Tomoaki Nakazono, Seiya Ishino, Yoshio Sakurai, Takashi Kitsukawa, Jong Won Lee, Hyunjung Lee, Min Whan Jung, Cecilia Babul, Pedro E. Maldonado, Kazutaka Takahashi, Fritzie I. Arce-McShane, Callum F. Ross, Barry J. Sessle, Nicholas G. Hatsopoulos, Thomas Brochier, Alexa Riehle, Paul Chorley, Sonja Grün, Hisao Nishijo, Satoe Ichihara-Takeda, Shintaro Funahashi, Keisetsu Shima, Hajime Mushiake, Yukako Yamane, Hiroshi Tamura, Ichiro Fujita, Naoko Inaba, Kenji Kawano, Sergei Kurkin, Kikuro Fukushima, Kiyoshi Kurata, Masato Taira, Ken Ichiro Tsutsui, Tadashi Ogawa, Hidehiko Komatsu, Kowa Koida, Keisuke Toyama, Barry J. Richmond, and Shigeru Shinomoto.  Similarity in Neuronal Firing Regimes across Mammalian Species.  Journal of Neuroscience, 36(21):5736–5747, may 2016.
  * [44] Nicholas A. Steinmetz, Cagatay Aydin, Anna Lebedeva, Michael Okun, Marius Pachitariu, Marius Bauza, Maxime Beau, Jai Bhagat, Claudia Böhm, Martijn Broux, Susu Chen, Jennifer Colonell, Richard J. Gardner, Bill Karsh, Fabian Kloosterman, Dimitar Kostadinov, Carolina Mora-Lopez, John O’Callaghan, Junchol Park, Jan Putzeys, Britton Sauerbrei, Rik J.J. van Daal, Abraham Z. Vollan, Shiwei Wang, Marleen Welkenhuysen, Zhiwen Ye, Joshua T. Dudman, Barundeb Dutta, Adam W. Hantman, Kenneth D. Harris, Albert K. Lee, Edvard I. Moser, John O’Keefe, Alfonso Renart, Karel Svoboda, Michael Häusser, Sebastian Haesler, Matteo Carandini, and Timothy D. Harris.  Neuropixels 2.0: A miniaturized high-density probe for stable, long-term brain recordings.  Science, 372(6539), apr 2021.
  * [45] Michael E. Rule, Timothy O’Leary, and Christopher D. Harvey.  Causes and consequences of representational drift.  Current Opinion in Neurobiology, 58:141–147, oct 2019.
  * [46] Joan del Castillo and Marta Pérez-Casany.  Overdispersed and underdispersed Poisson generalizations.  Journal of Statistical Planning and Inference, 134(2):486–500, oct 2005.
  * [47] Ian H. Stevenson.  Omitted Variable Bias in GLMs of Neural Spiking Activity.  Neural Computation, 30(12):3227–3258, oct 2018.
  * [48] Robbe L.T. Goris, J. Anthony Movshon, and Eero P. Simoncelli.  Partitioning neuronal variability.  Nature Neuroscience 2014 17:6, 17(6):858–865, apr 2014.
  * [49] V. J. Uzzell and E. J. Chichilnisky.  Precision of spike trains in primate retinal ganglion cells.  Journal of Neurophysiology, 92(2):780–789, aug 2004.
  * [50] André A. Fenton and Robert U. Muller.  Place cell discharge is extremely variable during individual passes of the rat through the firing field.  Proceedings of the National Academy of Sciences of the United States of America, 95(6):3182–3187, 1998.
  * [51] Riccardo Barbieri, Michael C Quirk, Loren M Frank, Matthew A Wilson, and Emery N Brown.  Construction and analysis of non-Poisson stimulus-response models of neural spiking activity.  Journal of Neuroscience Methods, 105(1):25–37, jan 2001.
  * [52] Abed Ghanbari, Christopher M. Lee, Heather L. Read, and Ian H. Stevenson.  Modeling stimulus-dependent variability improves decoding of population neural responses.  Journal of Neural Engineering, 16(6), 2019.

##  Appendix A Quantifying Uncertainties

After convergence, we have an approximation of the log-posterior
$P(\bm{\theta}_{t}|\bm{Y})\approx
N(\bm{\theta}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})$, and we can use this
approximation to quantify the uncertainty about the CMP parameters, as well as
about the mean rate at each time.

The CMP parameters are log-normal distributed. Let
$\bm{Z}_{it}=\begin{pmatrix}\bm{x}^{\prime}_{it}&\bm{0}\\\
\bm{0}&\bm{g}^{\prime}_{it}\end{pmatrix}$, then
$(\lambda_{it},\nu_{it})^{\prime}=\exp(\bm{Z}_{it}\bm{\theta}_{t})\sim\text{Lognormal}_{2}(\bm{Z}_{it}\bm{\mu}_{t},\bm{Z}_{it}\bm{\Sigma}_{t}\bm{Z}^{\prime}_{it})$.
Denote the variance of CMP parameters as $\bm{V}_{it}$., where for
$\bm{a}=\bm{Z\mu}$ and $\bm{S}=\bm{Z\Sigma Z^{\prime}}$,
$[V]_{mn}=e^{a_{m}+a_{n}+\frac{1}{2}(S_{mm}+S_{nn})}(e^{S_{mn}}-1)$

The conditional mean firing rate is $\delta_{it}=E(Y_{it})$, whose variance
can be calculated by the Delta method:

$\displaystyle\widehat{Var}(\delta_{it})$
$\displaystyle=\begin{pmatrix}\frac{\partial\delta_{it}}{\partial\lambda_{it}}&\frac{\partial\delta_{it}}{\partial\nu_{it}}\end{pmatrix}\bm{V}_{it}\begin{pmatrix}\frac{\partial\delta_{it}}{\partial\lambda_{it}}\\\
\frac{\partial\delta_{it}}{\partial\nu_{it}}\end{pmatrix}$ (6)
$\displaystyle\frac{\partial\delta_{it}}{\partial\lambda_{it}}$
$\displaystyle=\frac{\partial^{2}\log
Z_{it}}{\partial\log\lambda_{it}\partial\lambda_{it}}=\frac{Var(Y_{it})}{\lambda_{it}}$
(7) $\displaystyle\frac{\partial\delta_{it}}{\partial\nu_{it}}$
$\displaystyle=\frac{\partial^{2}\log
Z_{it}}{\partial\log\lambda_{it}\partial\nu_{it}}=-Cov(Y_{it},\log Y_{it}!)$
(8)

We can calculate the moments as in Appendix B, or we can use simpler
approximations $E(Y)=\lambda^{1/\nu}-\frac{\nu-1}{2\nu}$ when $\nu\leq 1$ or
$\lambda>10^{\nu}$. Then
$\frac{\partial\delta_{it}}{\partial\lambda_{it}}\approx\frac{1}{\nu_{it}}\lambda_{it}^{1/\nu_{it}-1}$
and
$\frac{\delta_{it}}{\nu_{it}}\approx-\frac{\lambda_{it}^{1/\nu_{it}}\log\lambda_{it}}{\nu_{it}^{2}}-\frac{1}{2\nu_{it}^{2}}$.

##  Appendix B Moments approximation for Conway-Maxwell Poisson distribution

To estimate the state-vector for the dynamic CMP model, we need to find first
and second moments for $Y$ and $\log Y!$. For $Y\sim CMP(\lambda,\nu)$,

$\displaystyle Z(\lambda,\nu)$
$\displaystyle=\sum_{k=0}^{\infty}\frac{\lambda^{k}}{(k!)^{\nu}}$ (9)
$\displaystyle E(Y)$ $\displaystyle=\frac{\partial\log
Z}{\partial\log\lambda}=\frac{1}{Z}\sum_{k=0}^{\infty}\frac{k\lambda^{k}}{(k!)^{\nu}}$
$\displaystyle Var(Y)$ $\displaystyle=\frac{\partial^{2}\log
Z}{\partial(\log\lambda)^{2}}=\frac{1}{Z}\sum_{k=0}^{\infty}\frac{k^{2}\lambda^{k}}{(k!)^{\nu}}-E^{2}(Y)$
$\displaystyle E(\log Y!)$ $\displaystyle=-\frac{\partial\log
Z}{\partial\nu}=\frac{1}{Z}\sum_{k=0}^{\infty}\frac{(\log
k!)\lambda^{k}}{(k!)^{\nu}}$ $\displaystyle Var(\log Y!)$
$\displaystyle=\frac{\partial^{2}\log
Z}{\partial\nu^{2}}=\frac{1}{Z}\sum_{k=0}^{\infty}\frac{(\log
k!)^{2}\lambda^{k}}{(k!)^{\nu}}-E^{2}(\log Y!)$ $\displaystyle Cov(Y,\log Y!)$
$\displaystyle=-\frac{\partial^{2}\log
Z}{\partial\log\lambda\partial\nu}=\frac{1}{Z}\sum_{k=0}^{\infty}\frac{(\log
k!)k\lambda^{k}}{(k!)^{\nu}}-E(\log Y!)E(Y)$

Generally, these moments can be calculated by truncated summation.

However, when $\lambda\geq 2$ and $\nu\leq 1$, we need many steps for accurate
approximation. In this case, we make use of a previous asymptotic results
([19, 24]) for efficient calculation. Let
$\alpha=\lambda^{1/\nu},c_{1}=\frac{\nu^{2}-1}{24}$ and
$c_{2}=\frac{\nu^{2}-1}{48}+\frac{c_{1}^{2}}{2}$,

$Z(\lambda,\nu)=\frac{e^{\nu\alpha}}{\lambda^{\frac{\nu-1}{2\nu}}(2\pi)^{\frac{\nu-1}{2}}\sqrt{\nu}}(1+c_{1}(\nu\alpha)^{-1}+c_{2}(\nu\alpha)^{-2}+\mathcal{O}(\lambda^{-3/\nu}))$
(10)

Then the moments can be calculated as follows:

$\displaystyle E(Y)$
$\displaystyle=\alpha-\frac{\nu-1}{2\nu}-\frac{\nu^{2}-1}{24\nu^{2}}\alpha^{-1}-\frac{\nu^{2}-1}{24\nu^{3}}\alpha^{-2}+\mathcal{O}(\alpha^{-3})$
(11) $\displaystyle Var(Y)$
$\displaystyle=\frac{\alpha}{\nu}+\frac{\nu^{2}-1}{24\nu^{3}}\alpha^{-1}+\frac{\nu^{2}-1}{12\nu^{4}}\aleph^{-2}+\mathcal{O}(\alpha^{-3})$
$\displaystyle E(\log Y!)$
$\displaystyle=\alpha\left(\frac{\log\lambda}{\nu}-1\right)+\frac{\log\lambda}{2\nu^{2}}+\frac{1}{2\nu}+\frac{\log
2\pi}{2}$
$\displaystyle-\frac{\alpha^{-1}}{24}\left(1+\frac{1}{\nu^{2}}+\frac{\log\lambda}{\nu}-\frac{\log\lambda}{\nu^{3}}\right)$
$\displaystyle-\frac{\alpha^{-2}}{24}\left(\frac{1}{\nu^{3}}+\frac{\log\lambda}{\nu^{2}}-\frac{\log\lambda}{\nu^{4}}\right)+\mathcal{O}(\alpha^{-3})$
$\displaystyle Var(\log Y!)$
$\displaystyle=\frac{\alpha(\log\lambda)^{2}}{\nu^{3}}+\frac{\log\lambda}{\nu^{3}}+\frac{1}{2\nu^{2}}$
$\displaystyle+\frac{\alpha^{-1}}{24\nu^{5}}[-2\nu^{2}+4\nu\log\lambda+(-1+\nu^{2})(\log\lambda)^{2}]$
$\displaystyle+\frac{\alpha^{-2}}{24\nu^{6}}[-3\nu^{2}-2\nu(-3+\nu^{2})\log\lambda+2(-1+\nu^{2})(\log\lambda)^{2}]+\mathcal{O}(\alpha^{-3})$
$\displaystyle Cov(Y,\log Y!)$
$\displaystyle=\frac{\alpha\log\lambda}{\nu^{2}}+\frac{1}{2\nu^{2}}+\frac{\alpha^{-1}}{24}\left(\frac{2}{\nu^{3}}+\frac{\log\lambda}{\nu^{2}}-\frac{\log\lambda}{\nu^{4}}\right)$
$\displaystyle-\frac{1}{24\alpha^{2}}\left(\frac{1}{\nu^{2}}-\frac{3}{\nu^{4}}-\frac{2\log\lambda}{\nu^{3}}+\frac{2\log\lambda}{\nu^{5}}\right)+\mathcal{O}(\alpha^{-3})$

##  Appendix C Gradient and Hessian of the log-posterior

We estimate the state vector by maximizing the log-posterior with Newton-
Raphson updates. Denote $f=P(\bm{\Theta|Y})$, the $(k+1)$-th update of NR
algorithm is
$\bm{\Theta}^{(k+1)}=\bm{\Theta}^{(k)}+[\nabla\nabla_{\bm{\Theta}^{(k)}}f]^{-1}\nabla_{\bm{\Theta}^{(k)}}f$.
The gradient is:

$\displaystyle\nabla_{\bm{\Theta}}f$ $\displaystyle=\left[\left(\frac{\partial
f}{\partial\bm{\theta}_{1}}\right)^{\prime},\ldots,\left(\frac{\partial
f}{\partial\bm{\theta}_{T}}\right)^{\prime}\right]^{\prime}$ (12)
$\displaystyle\frac{\partial f}{\partial\bm{\theta}_{1}}$
$\displaystyle=\frac{\partial
l_{1}}{\partial\bm{\theta}_{1}}-\bm{Q}_{0}^{-1}(\bm{\theta}_{1}-\bm{\theta}_{0})+\bm{F^{\prime}Q}^{-1}(\bm{\theta}_{2}-\bm{F\theta}_{1})$
$\displaystyle\frac{\partial f}{\partial\bm{\theta}_{t}}$
$\displaystyle=\frac{\partial
l_{t}}{\partial\bm{\theta}_{t}}-\bm{Q}^{-1}(\bm{\theta}_{t}-\bm{F\theta}_{t-1})+\bm{F^{\prime}Q}^{-1}(\bm{\theta}_{t+1}-\bm{F\theta}_{t})$
$\displaystyle\frac{\partial f}{\partial\bm{\theta}_{T}}$
$\displaystyle=\frac{\partial
l_{T}}{\partial\bm{\theta}_{T}}-\bm{Q}^{-1}(\bm{\theta}_{T}-\bm{F\theta}_{T-1})$
$\displaystyle\frac{\partial l_{t}}{\partial\bm{\theta}_{t}}$
$\displaystyle=\sum_{i=1}^{n}\begin{pmatrix}{\left(y_{it}-E\left(Y_{it}\right)\right)\bm{x}}_{it}\\\
{\nu_{it}\left(E\left(\log{Y_{it}!}\right)-\log{y_{it}!}\right)\bm{g}}_{it}\\\
\end{pmatrix}$

The Hessian is:

$\displaystyle\nabla\nabla_{\bm{\Theta}}f$
$\displaystyle=\begin{pmatrix}\frac{\partial^{2}f}{\partial\bm{\theta}_{1}\partial\bm{\theta}_{1}^{{}^{\prime}}}&\bm{F}^{\prime}\bm{Q}^{-1}&\bm{0}&\cdots&\bm{0}\\\
\bm{Q}^{-1}\bm{F}&\frac{\partial^{2}f}{\partial\bm{\theta}_{2}\partial\bm{\theta}_{2}^{{}^{\prime}}}&\bm{F}^{\prime}\bm{Q}^{-1}&\cdots&\vdots\\\
\bm{0}&\bm{Q}^{-1}\bm{F}&\frac{\partial^{2}f}{\partial\bm{\theta}_{3}\partial\bm{\theta}_{3}^{{}^{\prime}}}&\cdots&\vdots\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\bm{0}&\cdots&\cdots&\cdots&\frac{\partial^{2}f}{\partial\bm{\theta}_{T}\partial\bm{\theta}_{T}^{{}^{\prime}}}\end{pmatrix}$
(13)
$\displaystyle\frac{\partial^{2}f}{\partial\bm{\theta}_{1}\partial\bm{\theta}_{1}^{{}^{\prime}}}$
$\displaystyle=\frac{\partial^{2}l_{1}}{\partial\bm{\theta}_{1}\partial\bm{\theta}_{1}^{{}^{\prime}}}-\bm{Q}_{0}^{-1}-\bm{F}^{\prime}\bm{Q}^{-1}\bm{F}$
$\displaystyle\frac{\partial^{2}f}{\partial\bm{\theta}_{t}\partial\bm{\theta}_{t}^{{}^{\prime}}}$
$\displaystyle=\frac{\partial^{2}l_{t}}{\partial\bm{\theta}_{t}\partial\bm{\theta}_{t}^{{}^{\prime}}}-\bm{Q}^{-1}-\bm{F}^{\prime}\bm{Q}^{-1}\bm{F}$
$\displaystyle\frac{\partial^{2}f}{\partial\bm{\theta}_{T}\partial\bm{\theta}_{T}^{{}^{\prime}}}$
$\displaystyle=\frac{\partial^{2}l_{T}}{\partial\bm{\theta}_{T}\partial\bm{\theta}_{T}^{{}^{\prime}}}-\bm{Q}^{-1}$

, where

$\displaystyle\frac{\partial^{2}l_{t}}{\partial\bm{\theta}_{t}\partial\bm{\theta}_{t}^{{}^{\prime}}}$
$\displaystyle=\sum_{i=1}^{n}\begin{pmatrix}\bm{A}_{it}&\bm{B}_{it}\\\
\bm{B}^{\prime}_{it}&\bm{C}_{it}\end{pmatrix}$ (14) $\displaystyle\bm{A}_{it}$
$\displaystyle=Var(Y_{it})\bm{x}_{it}\bm{x}^{\prime}_{it}$
$\displaystyle\bm{B}_{it}$ $\displaystyle=-\nu_{it}Cov(Y_{it},\log
Y_{it}!)\bm{x}_{it}\bm{g}^{\prime}_{it}$ $\displaystyle\bm{C}_{it}$
$\displaystyle=\nu_{it}[\nu_{it}Var(\log Y_{it})-E(\log Y_{it}!)+\log
y_{it}!]\bm{g}_{it}\bm{g}^{\prime}_{it}$

When $\log y_{it}!\ll E(\log Y_{it}!)$, the Hessian may be ill-conditioned or
even positive-definite. To ensure the robustness, do Fisher scoring, i.e.
replace the observed information $-\nabla\nabla_{\bm{\Theta}}f$ by the
expected information $E(-\nabla\nabla_{\bm{\Theta}}f)$, so that
$\bm{C}_{it}=\nu_{it}^{2}Var(\log Y_{it}!)\bm{g}_{it}\bm{g}^{\prime}_{it}$