text:
\section{Introduction}

A key task facing natural or artificial agents acting in the real
world is that of causally inferring a hidden dynamic state based on
partial noisy observations. This task, referred to as \emph{filtering}
in the engineering literature, has been extensively studied since
the 1960s, both theoretically and practically (e.g., \cite{AndMoo05}).
For the linear setting with Gaussian noise and quadratic cost, the
solution is well known both for discrete and continuous times, leading
to the celebrated Kalman and the Kalman-Bucy filters \cite{KalBuc61,Kalman60},
respectively. In these cases the exact posterior state distribution
is Gaussian, resulting in closed form recursive update equations for
the mean and variance, yielding finite-dimensional filters. However,
beyond some very specific settings \cite{Daum05}, the optimal filter
is infinite-dimensional and impossible to compute in closed form,
requiring either approximate analytic techniques (e.g., the extended
Kalman filter (e.g., \cite{AndMoo05}), the unscented filter \cite{JulUhl00},
the cubature filter \cite{Arasaratnam2009}) or numerical procedures
(e.g., particle filters \cite{DouJoh09}). The latter usually require
time discretization and a finite number of particles, resulting in
loss of precision in a continuous time context.

The present work is motivated by the increasingly available data from
neuroscience, where the spiking activity of sensory neurons gives
rise to Point Process (PP)-like activity, which must be further analyzed
and manipulated by networks of neurons, in order to estimate attributes
of the external environment (e.g., the location and velocity of an
object), or to control the body, as in motor control, based on visual
and proprioceptive sensory inputs. In both these cases the system
faces the difficulty of assessing, in real-time, the environmental
state through a large number of simple, restricted and noisy sensors
(e.g., \cite{DayAbb05}). Example of such sensory cells are retinal
ganglion cells, auditory (cochlear) cells and proprioceptive stretch
receptors. In all these cases, information is conveyed to higher brain
areas through sequences of sharp pulses (spikes) delivered by the
responses of large numbers of sensory cells (about a million such
cells in the visual case). Each sensory cell is usually responsive
to a narrow set of attributes of the external stimuli (e.g., positions
in space, colors, frequencies etc.). Such cells are characterized
by \emph{tuning functions} or \emph{tuning curves}, with differential
responses to attributes of the external stimulus. For example, a visual
cell could respond with maximal probability to a stimulus at a specific
location in space and with reduced probability to stimuli distanced
from this point. An auditory cell could respond strongly to a certain
frequency range and with diminished responses to other frequencies,
etc. In all these cases, the actual firing of cells is random (due
to stochastic elements in the neurons, e.g., ion channels and synapses),
and can be described by PP with rate determined by the input and by
the cell's tuning function \cite{DayAbb05}. A particularly important
and ubiquitous phenomenon in biological sensory systems is \emph{sensory
adaptation}, which is the stimulus-dependent modification of system
parameters in a direction enhancing performance. Such changes usually
take place through the modification of tuning function properties,
e.g., the narrowing of tuning curve widths \cite{Benucci2009,GutZheKnu02},
the change of tuning curve heights \cite{Benucci2013}, or the shift
of the center position of tuning curves \cite{HarMcAlp04}. In the
sequel we refer to the process of setting the neurons' sensory tuning
functions as \emph{encoding}, and to the process of reconstructing
the state based on the PP observations as \emph{decoding}.

Inferring the hidden state under such circumstances has been widely
studied within the Computational Neuroscience literature, mostly for
static stimuli, homogeneous and equally-spaced tuning functions, and
using various approximations to the reconstruction error, such as
the Fisher information. In this work we are interested in setting
up a general framework for PP filtering in continuous time, and establishing
closed-form analytic expressions for an approximately optimal filter
(see \cite{BobMeiEld09,SusMeiOpp11,SusMeiOpp13} for previous work
in related, but more restrictive, settings). We aim to characterize
the nature of near-optimal encoders, namely to determine the structure
of the tuning functions for optimal state inference. A significant
advantage of the closed form expressions over purely numerical techniques
is the insight and intuition that is gained from them about \textit{qualitative}
aspects of the system. Moreover, the leverage gained by the analytic
computation contributes to reducing the variance inherent to Monte
Carlo approaches. Note that while this work has been motivated by
neuroscience, it should be viewed as a contribution to the theory
of point process filtering.

Technically, given the intractable infinite-dimensional nature of
the posterior distribution, we use a projection method replacing the
full posterior at each point in time by a projection onto a simple
family of distributions (Gaussian in our case). This approach, originally
developed in the Filtering literature \cite{Maybeck79,BriHanLeg98},
and termed Assumed Density Filtering (ADF), has been successfully
used more recently in Machine Learning \cite{Opper98,Minka01}. We
are aware of a single previous work using ADF in the context of point
process observations (\cite{SusemihlThesis2014}), where it was used
to optimize encoding by a population of two neurons.

Filtering PP observations based on multi-variate dynamic states in
continuous time has received far less attention in the engineering
literature than filtering based on more standard observations. While
a stochastic PDE for the infinite-dimensional posterior state distribution
can be derived under general conditions \cite{Snyder1972} (see also
\cite{Segall1976}), these equations are intractable in general, and
cannot even be qualitatively analyzed. This led \cite{RhoSny1977}
to consider a special case where the rate functions for the PP were
homogeneous Gaussians, leading to a posterior Gaussian distribution
of the state, giving rise to simple exact filtering equations for
the posterior mean and variance of the finite-dimensional posterior
state distribution. However, in our present setting we are motivated
to study inhomogeneous rate functions which lead to non-Gaussian infinite-dimensional
posterior distributions necessitating the introduction of finite-dimensional
approximations. Beyond neuroscience, PP based filtering has been used
for position sensing and tracking in optical communication \cite[sec. 4]{Snyder1977},
control of computer communication networks \cite{Segall1978}, queuing
\cite{Bremaud81} and econometrics \cite{Frey2001}, although our
main motivation for studying optimal sensory encoding arises from
neuroscience. Based on this motivation, we sometimes refer to point
events as ``spikes''.

The main contributions of the paper are the following: \emph{(i)}
Derivation of closed form recursive expressions for the continuous
time posterior mean and variance within the ADF approximation, allowing
for the incorporation of distributional assumptions over sensory variables
(going beyond homogeneity assumptions used so far). \emph{(ii)} Demonstrating
the quality of the ADF approximation by comparison to state of the
art particle filtering simulation methods. \emph{(iii)} Characterization
of optimal adaptation (encoding) for sensory cells in a more general
setting than hitherto considered. \emph{(iv)} Demonstrating the interesting
interplay between prior information and neuronal firing, showing how
in certain situations, the absence of spikes can be highly informative.
Preliminary results discussed in this paper were presented at a conference
\cite{Harel2015}. The present paper provides a rigorous formulation
of the mathematical framework and closed-form expressions for cases
that were not discussed in \cite{Harel2015}, including finite, possibly
heterogeneous, mixtures of Gaussian and uniform population.

\section{Problem Formulation}

\subsection{Heuristic formulation}

\begin{figure}
\bigskip{}

\includegraphics[width=1\columnwidth]{setting}

\bigskip{}
\caption{Problem setting}
\label{setting}
\end{figure}

We consider a dynamical system with state $X_{t}\in\mathbb{R}^{n}$,
for $t\in[0,\infty)$, observed through an array of sensors. Heuristically,
we assume the $i$-th sensor generates a random point in response
to $X_{t}=x$ with probability $\lambda_{t}\left(x;y_{i}\right)dt$
in the time interval $[t,t+dt)$, where $y_{i}\in\mathbf{Y}$ is a
parameter of the sensor. For example, in the case of auditory neurons,
$y_{i}$ may be the tuning curve center of the $i$-th neuron, corresponding
to the frequency for which it responds with the highest firing rate.
In this case the space $\mathbf{Y}$ would be the space of frequencies.

We assume that given $X_{t}=x$, points are generated independently
of the past of $X$ and of previously generated points from all sensors.
Each generated point is ``marked'' with the parameter $y_{i}$ of
the sensor that generated it. We therefore describe the observations
by a marked point process, i.e., a random counting measure $N$ on
$[0,\infty)\times\mathcal{\mathbf{Y}}$ where $\mathbf{Y}$ is the
\emph{mark space}, and the observations available at time $t$ are
the restriction of $N$ to the set $[0,t]\times\mathcal{\mathbf{Y}}$.
Such a process may also be described as a vector of (unmarked) point
processes, one for each sensor; however, we adopt the marked point
process view to allow taking the limit of infinitely many sensors,
as described below. Figure \ref{setting} illustrates this setting
in a biological context, where the sensors are neurons and the points
events are spikes (action potentials).

We denote by $f$ the counting measure of sensor marks over the mark
space $\mathbf{Y}$, i.e., $f=\sum_{i}\delta_{y_{i}}$, where $\delta_{y}$
is the Dirac delta measure at $y$. Then the marked process $N_{t}$
has the measure-valued random process $\lambda_{t}\left(X_{t};y\right)f\left(dy\right)$
as its \emph{intensity kernel} (see, e.g., \cite{Bremaud81}, Chapter
VIII), meaning that the rate (or \emph{intensity}) of points with
marks in a set $Y\subseteq\mathbf{Y}$ at time $t$ is
\[
\sum_{i:y_{i}\in Y}\lambda_{t}\left(X_{t};y_{i}\right)=\int_{Y}\lambda_{t}\left(X_{t};y\right)f\left(dy\right).
\]
Thus, the dynamics of $N$ may be described heuristically by means
of the intensity kernel as
\begin{equation}
\E\left[N\left(dt,dy\right)|X_{\left[0,t\right]},N_{\left[0,t\right]}\right]=\lambda_{t}\left(X_{t};y\right)f\left(dy\right)dt.\label{eq:rate-heuristic}
\end{equation}
Continuing the auditory example, the expected instantaneous total
rate of spikes from all neurons which are tuned to a frequency within
the range $\left[y_{\mathrm{min}},y_{\mathrm{max}}\right]$ (conditioned
on the history of the external state and of neural firing) would be
given by the integral of $\lambda_{t}\left(X_{t};y\right)f\left(dy\right)$
over this range of frequencies.

We generalize this model by allowing $f$ to be an arbitrary measure
on $\mathbf{Y}$ (not necessarily discrete). A continuous measure
may be useful as an approximate model for the response of a large
number of sensors, which may be easier to analyze than its discrete
counterpart. For example, a simplifying assumption in some previous
works is that tuning curve centers are located on an equally-spaced
grid (e.g., \cite{RhoSny1977,YaeMei10,Susemihl2014}). In the limit
of infinitely many neurons, this is equivalent to taking $f$ to be
the Lebesgue measure\textemdash the uniform ``distribution'' over
the entire space.

\subsection{Rigorous formulation}

We now proceed to describe our model more rigorously. We assume the
state $X_{t}\in\mathbb{R}^{n}$ obeys a Stochastic Differential Equation
(SDE)
\begin{equation}
dX_{t}=A\left(X_{t}\right)dt+B\left(U_{t}\right)dt+D\left(X_{t}\right)dW_{t},\quad\left(t\geq0\right),\label{eq:dynamics}
\end{equation}
where $A\left(\cdot\right),B\left(\cdot\right),D\left(\cdot\right)$
are arbitrary functions such that the SDE has a unique solution, $U$
is a control process and $W$ is a Wiener process. The initial condition
$X_{0}$ is assumed to have a continuous distribution with a known
density.

The observation process $N$ is a marked point process with marks
from a measurable space $\left(\mathbf{Y},\mathcal{Y}\right)$, i.e.,
a random counting measure on $[0,\infty)\times\mathbf{Y}$. We use
the notation $N_{t}\left(Y\right)\triangleq N\left([0,t]\times Y\right)$
for $Y\in\mathcal{Y}$, and when $Y=\mathbf{Y}$ we omit it and write
simply $N_{t}$.

Denote by $\left(\Omega,\mathcal{F},P\right)$ the underlying probability
space, and by $\mathcal{X},\mathcal{N}$ the filtrations generated
by $X,N$ respectively, that is,
\begin{align*}
\mathcal{X}_{t} & \triangleq\sigma\left(X_{\left[0,t\right]}\right)=\sigma\left\{ X_{s}:s\leq t\right\} \\
\mathcal{N}_{t} & \triangleq\sigma\left(N_{\left[0,t\right]\times\mathbf{Y}}\right)=\sigma\left\{ N\left((a,b]\times Y\right):a<b\leq t,Y\in\mathcal{Y}\right\}
\end{align*}
 where $N_{\left[0,t\right]\times\mathbf{Y}}$ is the measure $N$
restricted to $\left[0,t\right]\times\mathbf{Y}$. We assume the control
process $U$ is adapted to $\mathcal{N}$. Let $\mathcal{B}_{t}=\mathcal{X}_{t}\vee\mathcal{N}_{t}$
be the smallest $\sigma$-algebra containing $\mathcal{X}_{t}$ and
$\mathcal{N}_{t}$. We assume that the process $N_{t}\left(Y\right)$
has an intensity kernel with respect to $\mathcal{B}$ given by $\lambda_{t}\left(X_{t};y\right)f\left(dy\right)$
for some known function $\lambda_{t}\left(x;y\right)$ and measure
$f\left(dy\right)$ on $\mathcal{Y}$, meaning that for each $Y\in\mathcal{Y}$,
$C_{t}^{Y}\triangleq\int_{0}^{t}\left(\int_{Y}\lambda_{s}\left(X_{s};y\right)f\left(dy\right)\right)ds$
is $\mathcal{B}_{t}$-predictable and $N_{t}\left(Y\right)-C_{t}^{Y}$
is a $\mathcal{B}_{t}$-martingale. This condition is the rigorous
equivalent of (\ref{eq:rate-heuristic}) (see, e.g., \cite{Segall1975-modelling}
for a discussion of this definition in the unmarked case). Note in
particular that the presence of feedback in (\ref{eq:dynamics}),
which may depend on the history of $N$, means that $N\left(Y\right)$
is not in general a doubly-stochastic Poisson process, which forces
us to resort to a more sophisticated definition. Our results are applicable
to this general setting, and not restricted in application to doubly-stochastic
Poisson processes.

We define
\begin{align*}
\hat{\lambda}_{t}\left(y\right) & \triangleq\E\left[\lambda_{t}\left(X_{t};y\right)|\mathcal{N}_{t}\right],\\
\hat{\lambda}_{t}^{f} & \triangleq\int\hat{\lambda}_{t}\left(y\right)f\left(dy\right).
\end{align*}
The measure $\hat{\lambda}_{t}\left(y\right)f\left(dy\right)$ is
the intensity kernel of $N$ with respect to its natural filtration
$\mathcal{N}$, i.e., $\int_{Y}\hat{\lambda}_{t}\left(y\right)f\left(dy\right)$
is the intensity of $N_{t}\left(Y\right)$ w.r.t. $\mathcal{N}$ for
each $Y\in\mathcal{Y}$ (see \cite{Segall1975-modelling}, Theorem
2), and $\hat{\lambda}_{t}^{f}$ is the intensity of the unmarked
process $N_{t}$. For $Y\in\mathcal{Y}$, the innovation process of
$N_{t}\left(Y\right)$ is therefore $N_{t}\left(Y\right)-\int_{0}^{t}\int_{Y}\hat{\lambda}_{t}\left(y\right)f\left(dy\right)ds$,
and accordingly we define the \emph{innovation measure} $I$ by
\begin{equation}
I\left(dt,dy\right)\triangleq N\left(dt,dy\right)-\hat{\lambda}_{t}\left(y\right)f\left(dy\right)dt.\label{eq:innovation}
\end{equation}

\subsection{Encoding and decoding\label{sub:Encoding-and-decoding}}

We consider the question of optimal encoding and decoding under the
above model. By \emph{decoding }we mean computing (exactly or approximately)
the full posterior distribution of $X_{t}$ given $\mathcal{N}_{t}$.
The problem of optimal encoding is then the problem of optimal sensor
configuration, i.e., finding the optimal rate function $\lambda_{t}\left(x;y\right)$
and population distribution $f\left(dy\right)$ so as to minimize
some performance criterion. We assume the $\lambda_{t}\left(x;y\right)$
and $f\left(dy\right)$ come from some parameterized family with parameter
$\phi$.

To quantify the performance of the encoding-decoding system, we summarize
the result of decoding using a single estimator $\hat{X}_{t}=\hat{X}_{t}\left(\mathcal{N}_{t}\right)$,
and define the Mean Square Error (MSE) as $\epsilon_{t}\triangleq\mathrm{trace}[(X_{t}-\hat{X}_{t})(X_{t}-\hat{X}_{t})^{T}]$.
We seek $\hat{X}_{t}$ and $\phi$ that solve $\min_{\phi}\min_{\hat{X}_{t}}\mathrm{E}\left[\epsilon_{t}\right]=\min_{\phi}\mathrm{E}[\min_{\hat{X}_{t}}\mathrm{E}[\epsilon_{t}|\mathcal{N}_{t}]]$.
The inner minimization problem in this equation is solved by the MSE-optimal
decoder, which is the posterior mean $\hat{X}_{t}=\mu_{t}\triangleq\mathrm{E}\left[X_{t}|\mathcal{N}_{t}\right]$.
The posterior mean may be computed from the full posterior obtained
by decoding. The outer minimization problem is solved by the optimal
encoder. Note that, although we assume a fixed parameter $\phi$ which
does not depend on time, the optimal value of $\phi$ for which the
minimum is obtained generally depends on the time $t$ where the error
is to be minimized. In principle, the encoding/decoding problem can
be solved for any value of $t$. In order to assess performance it
is convenient to consider the steady-state limit $t\to\infty$ for
the encoding problem.

Below, we approximately solve the decoding problem for any $t$, for
some specific forms of $\lambda_{t}$ and $f$. We then explore the
problem of choosing the steady-state optimal encoding parameters $\phi$
using Monte Carlo simulations in a biologically-motivated toy model.
Note that if decoding is exact, the problem of optimal encoding becomes
that of minimizing the expected posterior variance.

Having an efficient (closed-form) approximate filter allows performing
the Monte Carlo simulation at a significantly reduced computational
cost, relative to numerical methods such as particle filtering. The
computational cost is further reduced by averaging the computed posterior
variance across trials, rather than the squared error, thereby requiring
fewer trials. The mean of the posterior variance equals the MSE (of
the posterior mean), but has the advantage of being less noisy than
the squared error itself \textendash{} since by definition it is the
mean of the square error under conditioning on $\mathcal{N}_{t}$
.

\subsection{Special cases\label{sub:special-cases}}

To approximately solve the decoding problem in closed form, we focus
on the case of Gaussian sensors: each sensor is marked with $y=\left(h,\theta,R\right)$
where $h\in\mathbb{R}_{+},\theta\in\mathbb{R}^{m}$, $R\in\mathbb{R}^{m\times m}$
is positive semidefinite, $m\leq n$, and
\begin{equation}
\lambda_{t}\left(x;h,\theta,R\right)=h\exp\left(-\frac{1}{2}\left\Vert Hx-\theta\right\Vert _{R}^{2}\right),\label{eq:guassian-sensor}
\end{equation}
where $\left\Vert z\right\Vert _{M}^{2}\triangleq z^{T}Mz$, and $H\in\mathbb{R}^{m\times n}$
is a fixed matrix of full row rank, which maps the external state
from $\mathbb{R}^{n}$ to ``sensory coordinates'' in $\mathbb{R}^{m}$.
Here, $h$ is the tuning curve height, i.e., the maximum firing rate
of the sensor; $\theta$ is the tuning curve center, i.e., the stimulus
value in sensory coordinates which elicits the highest firing rate;
and $R$ is the precision matrix of the Gaussian response curve in
the sensory space $\mathbb{R}^{m}$. Following neuroscience terminology,
we refer to the tuning curve center $\theta$ as the sensor's \emph{preferred
stimulus. }The inclusion of the matrix $H$ allows using high-dimensional
models where only some dimensions are observed, for example when the
full state includes velocities but only locations are directly observable.
In the one-dimensional case, $R^{-1/2}$ is the tuning curve width.

We consider several special forms of the population distribution $f$$\left(dh,d\theta,dR\right)$
where we can bring the approximate filter to closed form:

\subsubsection{A single sensor}

$f\left(dy\right)=\delta_{y}\left(dy\right)$, where $y=\left(h,\theta,R\right)$.

\subsubsection{Uniform population}

Here $h,R$ are fixed across all sensors, and $f\left(dh',d\theta,dR'\right)=\delta_{h}\left(dh'\right)\delta_{R}\left(dR'\right)d\theta$.

\subsubsection{Gaussian population}

Abusing notation slightly, we write $f\left(dh',d\theta,dR'\right)=\delta_{h}\left(dh'\right)\delta_{R}\left(dR'\right)f\left(d\theta\right)$
where

\begin{align}
f\left(d\theta\right) & =\mathcal{N}\left(\theta;c,\Sigma_{\mathrm{pop}}\right)d\theta\nonumber \\
 & =\left(2\pi\right)^{-n/2}\left|\Sigma_{\mathrm{pop}}\right|^{-1/2}\exp\left(-\frac{1}{2}\left\Vert \theta-c\right\Vert _{\Sigma_{\mathrm{pop}}^{-1}}^{2}\right)d\theta,\label{eq:gaussian-f}
\end{align}
for fixed $c\in\mathbb{R}^{m}$, and positive definite $\Sigma_{\mathrm{pop}}$.

We take $f$ to be normalized, since any scaling of $f$ may be included
in the coefficient $h$ in (\ref{eq:guassian-sensor}), resulting
in the same point-process. Thus, when used to approximate a large
population of sensors, the coefficient $h$ would be proportional
to the number of sensors.

\subsubsection{Uniform population on an interval}

In this case we assume a scalar state, $n=m=1$, and

\begin{align}
f\left(d\theta\right) & =1_{\left[a,b\right]}\left(\theta\right)d\theta,\label{eq:interval}
\end{align}
where similarly to the Gaussian population case, $h$ and $R$ are
fixed. Unlike the Gaussian case, here we find it more convenient not
to normalize the distribution.

\subsubsection{Finite mixtures}

$f\left(dy\right)=\sum_{i}\alpha_{i}f_{i}\left(dy\right),$ where
each $f_{i}$ is of one of the above forms. This form is quite general:
it includes populations where $\theta$ is distributed according to
a Gaussian mixture, as well as heterogeneous populations with finitely
many different values of $R$\@. The resulting filter derived below
includes a term for each component of the mixture.

\section{Assumed Density Filtering}

\subsection{Exact filtering equations}

Let $P\left(\cdot,t\right)$ denote the posterior density of $X_{t}$
given $\mathcal{N}_{t}$, and $\mathrm{E}_{P}^{t}\left[\cdot\right]$
the posterior expectation given $\mathcal{N}_{t}$. The prior density
$P\left(\cdot,0\right)$ is assumed to be known.

The problem of filtering a diffusion process $X$ from a doubly stochastic
Poisson process driven by $X$ is formally solved in \cite{Snyder1972}.
The result is extended to general marked point processes in the presence
of feedback in \cite{RhoSny1977}, where the authors derive a stochastic
PDE for the posterior density, which in our setting takes the form
\begin{align}
dP\left(x,t\right)= & \mathcal{L}_{t}^{*}P\left(x,t\right)dt+P\left(x,t^{-}\right)\int_{y\in\mathbf{Y}}\frac{\lambda_{t^{-}}\left(x;y\right)-\hat{\lambda}_{t^{-}}\left(y\right)}{\hat{\lambda}_{t^{-}}\left(y\right)}I\left(dt,dy\right),\label{eq:rhodes-snyder}
\end{align}
where $\mathcal{L}_{t}$ is the state's posterior infinitesimal generator
(Kolmogorov's backward operator), defined as $\mathcal{L}_{t}f\left(x\right)=\lim_{\Delta t\to0^{+}}\left(\mathrm{E}\left[f\left(X_{t+\Delta t}\right)|\mathcal{B}_{t}\right]-f\left(x\right)\right)/\Delta t$,
$\mathcal{L}_{t}^{*}$ is $\mathcal{L}_{t}$'s adjoint operator (Kolmogorov's
forward operator), and $I$ is defined in (\ref{eq:innovation}).
Expressions including $t^{-}$ are to be interpreted as left limits\footnote{The formulation in \cite{RhoSny1977} does not include left limits,
since it uses conditioning only on the strict past, resulting in a
left-continuous posterior. Our definition of $\mc N_{t}$ includes
the current time $t$, following the convention used in \cite{Segall1975-modelling}
and others, so that left limits are required at spike times.}. Note that the discrete part of the measure $I\left(dt,dy\right)$,
namely $N\left(dt,dy\right)$, contributes a discontinuous change
in the posterior at each spike time. Also note that in this closed-loop
setting, the infinitesimal generator is itself a random operator,
due to its dependence on past observations through the control law.

The stochastic PDE (\ref{eq:rhodes-snyder}) is non-linear and non-local
(due to the dependence of $\hat{\lambda}_{t}\left(y\right)$ on $P\left(\cdot,t\right)$),
and therefore usually intractable. In \cite{RhoSny1977,Susemihl2014}
the authors consider linear dynamics with a Gaussian prior and Gaussian
sensors with centers distributed uniformly over the state space. In
this case, the posterior is Gaussian, and (\ref{eq:rhodes-snyder})
leads to closed-form ODEs for its moments. In our more general setting,
we can obtain exact equations for the posterior moments, as follows.

Let $\mu_{t}\triangleq\mathrm{E}_{P}^{t}X_{t},\tilde{X}_{t}\triangleq X_{t}-\mu_{t},\Sigma_{t}\triangleq\mathrm{E}_{P}^{t}[\tilde{X}_{t}\tilde{X}_{t}^{T}]$.
Using (\ref{eq:rhodes-snyder}), along with known results about the
form of the infinitesimal generator $\mathcal{L}_{t}$ for diffusion
processes (e.g. \cite{Oksendal2003}, Theorem 7.3.3), the first two
posterior moments can be shown to obey the following equations\footnote{see \cite{SusMeiOpp13} for derivation between spikes, or \cite{Segall1975-modelling}
for a derivation of (\ref{eq:mean}) via a different method}:
\begin{align}
d\mu_{t} & =\left(\mathrm{E}_{P}^{t}\left[A\left(X_{t}\right)\right]+B\left(U_{t}\right)\right)dt+\int_{\mathbf{Y}}\mathrm{E}_{P}^{t-}\left[\omega_{t^{-}}\left(y\right)X_{t^{-}}\right]I\left(dt,dy\right)\label{eq:mean}\\
d\Sigma_{t} & =\mathrm{E}_{P}^{t}\left[A\left(X_{t\left(y\right)}\right)\tilde{X}_{t}^{T}+\tilde{X}_{t}A\left(X_{t}\right)^{T}+D\left(X_{t}\right)D\left(X_{t}\right)^{T}\right]dt\nonumber \\
 & \quad+\int_{\mathbf{Y}}\mathrm{E}_{P}^{t^{-}}\left[\omega_{t^{-}}\left(y\right)\tilde{X}_{t^{-}}\tilde{X}_{t^{-}}^{T}\right]I\left(dt,dy\right)\nonumber \\
 & \quad-\int_{\mathbf{Y}}\mathrm{E}_{P}^{t^{-}}\left[\omega_{t^{-}}\left(y\right)X_{t^{-}}\right]\mathrm{E}_{P}^{t^{-}}\left[\omega_{t^{-}}\left(y\right)X_{t^{-}}^{T}\right]N\left(dt,dy\right)\label{eq:var}
\end{align}
where
\[
\omega_{t}\left(y\right)\triangleq\frac{\lambda_{t}\left(X_{t};y\right)}{\hat{\lambda}_{t}\left(y\right)}-1,
\]
and similarly we write $\omega_{t}\left(x;y\right)\triangleq\lambda_{t}\left(x;y\right)/\hat{\lambda}_{t}\left(y\right)-1$.

In contrast with the more familiar LQG problem and the Kalman-Bucy
filter, here the posterior variance is random, and is generally not
monotonically decreasing even when estimating a constant state. However,
noting that $\E\left[I\left(dt,dy\right)\right]=0$, we may observe
from (\ref{eq:var}) that for a constant state ($A=D=0$), the \emph{expected}
posterior variance $\E\left[\Sigma_{t}\right]$ is decreasing.

Equations (\ref{eq:mean})-(\ref{eq:var}) are written in terms of
the innovation process $I\left(dt,dy\right)$ and the original point
process. Although this formulation highlights the relation to the
Kalman-Bucy filter, we will find it useful to rewrite them in a different
form, as follows,
\begin{align*}
d\mu_{t} & =d\mu_{t}^{\pi}+d\mu_{t}^{\mathrm{c}}+d\mu_{t}^{N},\\
d\Sigma_{t} & =d\Sigma_{t}^{\pi}+d\Sigma_{t}^{\mathrm{c}}+d\Sigma_{t}^{N},
\end{align*}
where $d\mu_{t}^{\pi},d\Sigma_{t}^{\pi}$ are the \emph{prior terms},\emph{
}corresponding to $\mathcal{L}_{t}^{*}P\left(x,t\right)$ in (\ref{eq:rhodes-snyder}),
and the remaining terms are divided into \emph{continuous update terms}
$d\mu_{t}^{\mathrm{c}},d\Sigma_{t}^{\mathrm{c}}$ (multiplying $dt$)
and \emph{discontinuous update terms} $d\mu_{t}^{N},d\Sigma_{t}^{N}$
(multiplying $N\left(dt,dy\right)$). Using (\ref{eq:innovation}),
we find
\begin{align}
d\mu_{t}^{\pi} & =\mathrm{E}_{P}^{t}\left[A\left(X_{t}\right)\right]dt+B\left(U_{t}\right)dt\label{eq:mean-p}\\
d\Sigma_{t}^{\pi} & =\mathrm{E}_{P}^{t}\left[A\left(X_{t}\right)\tilde{X}_{t}^{T}+\tilde{X}_{t}A\left(X_{t}\right)^{T}+D\left(X_{t}\right)D\left(X_{t}\right)^{T}\right]dt\label{eq:var-p}
\end{align}
\begin{align}
d\mu_{t}^{\mathrm{c}} & =-\int_{\mathbf{Y}}\mathrm{E}_{P}^{t}\left[\omega_{t}\left(y\right)X_{t}\right]\hat{\lambda}_{t}\left(y\right)f\left(dy\right)dt\label{eq:mean-c}\\
d\Sigma_{t}^{\mathrm{c}} & =-\int_{\mathbf{Y}}\mathrm{E}_{P}^{t}\left[\omega_{t}\left(y\right)\tilde{X}_{t}\tilde{X}_{t}^{T}\right]\hat{\lambda}_{t}\left(y\right)f\left(dy\right)dt\label{eq:var-c}
\end{align}

\begin{align}
d\mu_{t}^{N} & =\int_{\mathbf{Y}}\mathrm{E}_{P}^{t-}\left[\omega_{t^{-}}\left(y\right)X_{t^{-}}\right]N\left(dt,dy\right)\label{eq:mean-N}\\
d\Sigma_{t}^{N} & =\int_{\mathbf{Y}}\bigg(\mathrm{E}_{P}^{t^{-}}\left[\omega_{t^{-}}\left(y\right)\tilde{X}_{t^{-}}\tilde{X}_{t^{-}}^{T}\right]\label{eq:var-N}\\
 & \quad-\mathrm{E}_{P}^{t^{-}}\left[\omega_{t^{-}}\left(y\right)X_{t^{-}}\right]\mathrm{E}_{P}^{t^{-}}\left[\omega_{t^{-}}\left(y\right)X_{t^{-}}^{T}\right]\bigg)N\left(dt,dy\right).\nonumber
\end{align}

The prior terms $d\mu_{t}^{\pi},d\Sigma_{t}^{\pi}$ represent the
known dynamics of $X$, and are the same terms appearing in the Kalman-Bucy
filter. These would be the only terms left if no measurements were
available, and would vanish for a static state. The continuous update
terms $d\mu_{t}^{\mathrm{c}},d\Sigma_{t}^{\mathrm{c}}$ represent
updates to the posterior between spikes that are not derived from
$X$'s dynamics, and therefore may be interpreted as corresponding
to information obtained from the absence of spikes. The discontinuous
update terms $d\mu_{t}^{N},d\Sigma_{t}^{N}$ contribute a change to
the posterior at spike times, depending on the spike's mark $y$,
and thus represent information obtained from the presence of a spike
as well as its associated mark.

Note that only the continuous update terms depend explicitly on the
population distribution $f$. Discontinuous updates depend on the
population distribution only indirectly through their influence on
the statistics of the point process $N$.

\subsection{ADF approximation}

While equations (\ref{eq:mean-p})-(\ref{eq:var-N}) are exact, they
are not practical, since they require computation of the full posterior
$\mathrm{E}_{P}^{t}\left[\cdot\right]$. To bring them to a closed
form, we use ADF with an assumed Gaussian density (see \cite{Opper98}
for details). Informally, this may be envisioned as integrating (\ref{eq:mean-p})-(\ref{eq:var-N})
while replacing the distribution $P$ by its approximating Gaussian
``at each time step''. The approximating Gaussian is obtained by
matching the first two moments of $P$ \cite{Opper98}. Note that
the solution of the resulting equations does not in general match
the first two moments of the exact solution, though it may approximate
it. Practically, the ADF approximation amounts to substituting the
normal distribution $\mathcal{N}(x;\mu_{t},\Sigma_{t})$ for $P(x,t)$
to compute the expectations in (\ref{eq:mean-p})-(\ref{eq:var-N}).

If the dynamics are linear, the prior updates (\ref{eq:mean-p})-(\ref{eq:var-p})
are easily computed in closed form after this substitution. Otherwise,
they may be approximated assuming the non-linear functions $A(x)$
and $D\left(x\right)D\left(x\right)^{T}$ may be written as power
series, using well-known results about the moments of Gaussian vectors.
The next sections are therefore devoted to the approximation of the
non-prior updates (\ref{eq:mean-c})-(\ref{eq:var-N}).

\subsection{Approximate form for Gaussian sensors}

We now proceed to apply the Gaussian ADF approximation $P\left(x,t\right)\approx\mathcal{N}\left(x;\mu_{t},\Sigma_{t}\right)$
to (\ref{eq:mean-c})-(\ref{eq:var-N}) in the case of Gaussian sensors
(\ref{eq:guassian-sensor}), deriving approximate filtering equations
written in terms of the population density $f\left(dh,d\theta,dR\right)$.
Abusing notation, from here on we use $\mu_{t},\Sigma_{t}$, and $P\left(x,t\right)$
to refer to the ADF approximation rather than to the exact values.

To evaluate the posterior of expectations in (\ref{eq:mean-c})-(\ref{eq:var-N})
we first simplify the expression
\begin{equation}
P\left(x,t\right)\omega_{t}\left(x;y\right)=\frac{P\left(x,t\right)\lambda_{t}\left(x;y\right)}{\int P\left(\xi,y\right)\lambda_{t}\left(\xi;y\right)d\xi}-P\left(x,t\right).\label{eq:P-omega}
\end{equation}
Using the Gaussian ADF approximation $P\left(x,t\right)=\mathcal{N}\left(x;\mu_{t},\Sigma_{t}\right)$
and (\ref{eq:guassian-sensor}), we find
\begin{align*}
P\left(x,t\right) & \lambda_{t}\left(x;h,\theta,R\right)=\\
 & =h\mathcal{N}\left(x;\mu_{t},\Sigma_{t}\right)\exp\left(-\frac{1}{2}\left\Vert Hx-\theta\right\Vert _{R}^{2}\right)\\
 & =\frac{h\exp\left(-\frac{1}{2}\left\Vert x-\mu_{t}\right\Vert _{\Sigma_{t}^{-1}}^{2}-\frac{1}{2}\left\Vert Hx-\theta\right\Vert _{R}^{2}\right)}{\sqrt{\left(2\pi\right)^{n}\left|\Sigma_{t}\right|}}\\
 & =\frac{h}{\sqrt{\left(2\pi\right)^{n}\left|\Sigma_{t}\right|}}\exp\left(-\frac{1}{2}\left\Vert H_{r}^{-1}\theta-\mu_{t}\right\Vert _{Q_{t}^{R}}^{2}-\frac{1}{2}\left\Vert x-\mu_{t}^{\theta}\right\Vert _{\Sigma_{t}^{-1}+H^{T}RH}^{2}\right),
\end{align*}
where $H_{r}^{-1}$ is any right inverse of $H$, and
\begin{eqnarray*}
Q_{t}^{R} & \triangleq & \Sigma_{t}^{-1}\left(\Sigma_{t}^{-1}+H^{T}RH\right)^{-1}H^{T}RH,\\
\mu_{t}^{\theta} & \triangleq & \left(\Sigma_{t}^{-1}+H^{T}RH\right)^{-1}\left(Q_{t}\mu_{t}+H^{T}R\theta\right).
\end{eqnarray*}
An application of the Woodbury identity establishes the relation $Q_{t}^{R}=H^{T}S_{t}^{R}H$,
where
\begin{equation}
S_{t}^{R}\triangleq\left(R^{-1}+H\Sigma_{t}H^{T}\right)^{-1},\label{eq:S}
\end{equation}
yielding
\begin{multline}
P\left(x,t\right)\lambda_{t}\left(x;h,\theta,R\right)\\
\quad=\frac{h}{\sqrt{\left(2\pi\right)^{n}\left|\Sigma_{t}\right|}}\exp\left(-\frac{1}{2}\left\Vert \theta-H\mu_{t}\right\Vert _{S_{t}^{R}}^{2}-\frac{1}{2}\left\Vert x-\mu_{t}^{\theta}\right\Vert _{\Sigma_{t}^{-1}+H^{T}RH}^{2}\right),\label{eq:post-times-rate}
\end{multline}
and by normalizing this Gaussian (see (\ref{eq:P-omega})) we find
that
\begin{align*}
P\left(x,t\right)\omega_{t}\left(x;h,\theta,R\right) & =\mathcal{N}\left(x;\mu_{t}^{\theta},\left(\Sigma_{t}^{-1}+H^{T}RH\right)^{-1}\right)-P\left(x,t\right)
\end{align*}
yielding
\begin{align*}
\mathrm{E}_{P}^{t}\left[\omega_{t}\left(h,\theta,R\right)X_{t}\right] & =\mu_{t}^{\theta}-\mu_{t}\\
\mathrm{E}_{P}^{t}\left[\omega_{t}\left(h,\theta,R\right)\tilde{X}_{t}\tilde{X}_{t}^{T}\right] & =\left(\Sigma_{t}^{-1}+H^{T}RH\right)^{-1}+\left(\mu_{t}-\mu_{t}^{\theta}\right)\left(\mu_{t}-\mu_{t}^{\theta}\right)^{T}-\Sigma_{t}.
\end{align*}
Substituting the definition of $\mu_{t}^{\theta}$ and simplifying
using the Woodbury identity yields
\begin{align*}
\mathrm{E}_{P}^{t} & \left[\omega_{t}\left(h,\theta,R\right)X_{t}\right]=-\Sigma_{t}H^{T}S_{t}^{R}\left(H\mu_{t}-\theta\right),\\
\mathrm{E}_{P}^{t} & \left[\omega_{t}\left(h,\theta,R\right)\tilde{X}_{t}\tilde{X}_{t}^{T}\right]=\Sigma_{t}H^{T}\left(S_{t}^{R}\left(H\mu_{t}-\theta\right)\left(H\mu_{t}-\theta\right)^{T}S_{t}^{R}-S_{t}^{R}\right)H\Sigma_{t}.
\end{align*}
Plugging this result into (\ref{eq:mean-c})-(\ref{eq:var-N}) yields

\begin{align}
d\mu_{t}^{\mathrm{c}} & =\Sigma_{t}H^{T}\int_{\mathbf{Y}}S_{t}^{R}\left(H\mu_{t}-\theta\right)\hat{\lambda}_{t}\left(h,\theta,R\right)f\left(dh,d\theta,dR\right)dt\label{eq:mean-c-gaussian}\\
d\Sigma_{t}^{\mathrm{c}} & =\Sigma_{t}H^{T}\int_{\mathbf{Y}}\left(S_{t}^{R}-S_{t}^{R}\left(H\mu_{t}-\theta\right)\left(H\mu_{t}-\theta\right)^{T}S_{t}^{R}\right)\nonumber \\
 & \qquad\qquad\qquad\times\hat{\lambda}_{t}\left(h,\theta,R\right)f\left(dh,d\theta,dR\right)H\Sigma_{t}dt\label{eq:var-c-gaussian}\\
d\mu_{t}^{N} & =\Sigma_{t^{-}}H^{T}\int_{\mathbf{Y}}S_{t^{-}}^{R}\left(\theta-H\mu_{t^{-}}\right)N\left(dt,d\left[h,\theta,R\right]\right)\label{eq:mean-N-guassian}\\
d\Sigma_{t}^{N} & =-\Sigma_{t^{-}}H^{T}\int_{\mathbf{Y}}S_{t^{-}}^{R}N\left(dt,d\left[h,\theta,R\right]\right)H\Sigma_{t^{-}}\label{eq:var-N-gaussian}
\end{align}

We also note that integrating (\ref{eq:post-times-rate}) over $x$
yields
\begin{align}
\hat{\lambda}_{t}\left(h,\theta,R\right) & =h\sqrt{\frac{\left|S_{t}^{R}\right|}{\left|R\right|}}\exp\left(-\frac{1}{2}\left\Vert \theta-H\mu_{t}\right\Vert _{S_{t}^{R}}^{2}\right),\label{eq:self-rate-gaussian}
\end{align}
where we computed the coefficient using the equality $\left|I+\Sigma_{t}H^{T}RH\right|=\left|R\right|/\left|S_{t}^{R}\right|$,
derived by application of Sylvester's determinant identity.

To gain some insight into these equations, consider the case $H=I$.
As seen from the discontinuous update equations (\ref{eq:mean-N-guassian})-(\ref{eq:var-N-gaussian}),
when a spike with mark $\left(h,\theta,R\right)$ occurs, the posterior
mean moves towards the location mark $\theta$, and the posterior
variance decreases (in the sense that $\Sigma_{t^{+}}-\Sigma_{t^{-}}$
is negative definite). Neither update depends on $h$.

In the scalar case $m=n=1$, with $H=1$, $\sigma_{t}^{2}=\Sigma_{t},\sigma_{\mathrm{r}}^{2}=R^{-1},a=A,d=D$,
the discontinuous update equations (\ref{eq:mean-N-guassian})-(\ref{eq:var-N-gaussian})
read
\begin{align}
d\mu_{t}^{N} & =\int_{\mathbf{Y}}\frac{\sigma_{t^{-}}^{2}}{\sigma_{t^{-}}^{2}+\sigma_{\mathrm{r}}^{2}}\left(\theta-\mu_{t^{-}}\right)N\left(dt,d\left[h,\theta,\sigma_{r}^{2}\right]\right)\label{eq:mean-N-guassian-1d}\\
d\Sigma_{t}^{N} & =-\int_{\mathbf{Y}}\frac{\sigma_{t^{-}}^{2}}{\sigma_{t^{-}}^{2}+\sigma_{\mathrm{r}}^{2}}\sigma_{t^{-}}^{2}N\left(dt,d\left[h,\theta,\sigma_{r}^{2}\right]\right)\label{eq:var-N-gaussian-1d}
\end{align}

The continuous mean update equation (\ref{eq:mean-c-gaussian}) also
admits an intuitive interpretation, in the case where all sensors
share the same precision matrix $R$, i.e. $f\left(d\left[h,\theta,R'\right]\right)=f\left(dh,d\theta\right)\delta_{R}\left(dR'\right)$.
In this case, the equation reads
\begin{align*}
d\mu_{t}^{\mathrm{c}} & =\Sigma_{t}H^{T}S_{t}^{R}\left(H\mu_{t}-\int_{\mathbf{Y}}\theta\nu_{t}\left(\theta\right)f\left(dh,d\theta\right)\right)\hat{\lambda}_{t}^{f}dt
\end{align*}
where $\nu_{t}\left(\theta\right)\triangleq\hat{\lambda}_{t}\left(h,\theta,R\right)/\hat{\lambda}_{t}^{f}$
(which does not depend on $h$). The normalized intensity kernel $\nu_{t}\left(\theta\right)f\left(dh,d\theta\right)$
may be interpreted heuristically as the distribution of the mark $\left(h,\theta\right)$
of the next event provided it occurs immediately. Thus, the absence
of events drives the posterior mean away from the expected location
of the next mark. The strength of this effect scales with $\hat{\lambda}_{t}^{f}$,
which is the rate of the unmarked process $N_{t}$ w.r.t.~its natural
filtration, i.e., the total expected rate of spikes given the firing
history. This behavior is qualitatively similar to the result obtained
in \cite{BobMeiEld09} for a finite population of sensors observing
a continuous-time finite-state Markov process, where the posterior
probability between spikes concentrates on states with lower total
firing rate.

\subsection{Closed form approximations in special cases}

Using (\ref{eq:self-rate-gaussian}), we now evaluate the continuous
update equations (\ref{eq:mean-c-gaussian})-(\ref{eq:var-c-gaussian})
for the specific forms of the population distribution $f\left(dy\right)$
listed in section \ref{sub:special-cases}. Note that the discontinuous
update equations (\ref{eq:mean-N-guassian})-(\ref{eq:var-N-gaussian})
do not depend on the population distribution $f$, and are already
in closed form.

\subsubsection{Single sensor}

The result for a single sensor with parameters $h,\theta,R$ is trivial
to obtain from (\ref{eq:mean-c-gaussian})-(\ref{eq:var-c-gaussian}),
yielding
\begin{align}
d\mu_{t}^{\mathrm{c}} & =\Sigma_{t}H^{T}S_{t}^{R}\left(H\mu_{t}-\theta\right)\hat{\lambda}_{t}^{f}dt,\label{eq:mean-c-single}\\
d\Sigma_{t}^{\mathrm{c}} & =\Sigma_{t^{-}}H^{T}\left(S_{t}^{R}-S_{t}^{R}\left(H\mu_{t}-\theta\right)\left(H\mu_{t}-\theta\right)^{T}S_{t}^{R}\right)\hat{\lambda}_{t}^{f}H\Sigma_{t^{-}}dt,\label{eq:var-c-single}
\end{align}
where $\hat{\lambda}_{t}^{f}=\hat{\lambda}_{t}\left(h,\theta,R\right)$
as given by (\ref{eq:self-rate-gaussian}), and $S_{t}^{R}$ is defined
in (\ref{eq:S}).

\subsubsection{Uniform population}

Here all sensors share the same height $h$ and precision $R$, whereas
the location parameter $\theta$ covers $\mathbb{R}^{m}$ uniformly,
i.e. $f\left(d\theta\right)=d\theta$. A straightforward calculation
from (\ref{eq:mean-c-gaussian})-(\ref{eq:var-c-gaussian}) and (\ref{eq:self-rate-gaussian})
yields
\begin{align}
d\mu_{t}^{\mathrm{c}} & =0,\quad d\Sigma_{t}^{\mathrm{c}}=0,\label{eq:mean-var-uniform}
\end{align}
in agreement with the (exact) result obtained in \cite{RhoSny1977},
where the filtering equations only include the prior term and the
discontinuous update term.

\subsubsection{Gaussian population}

Here $R,h$ are constant, and $f\left(d\theta\right)$ is given in
(\ref{eq:gaussian-f}). Using (\ref{eq:self-rate-gaussian}),
\begin{align*}
\hat{\lambda}_{t}\left(h,\theta,R\right)f\left(d\theta\right) & =h\sqrt{\frac{\left|S_{t}^{R}\right|}{\left|R\right|}}\exp\left(-\frac{1}{2}\left\Vert \theta-H\mu_{t}\right\Vert _{S_{t}^{R}}^{2}\right)\mathcal{N}\left(\theta;c,\Sigma_{\mathrm{pop}}\right)d\theta.
\end{align*}
An analogous computation to the derivation of (\ref{eq:post-times-rate})
and (\ref{eq:self-rate-gaussian}) above yields
\begin{align}
\hat{\lambda}_{t}\left(h,\theta,R\right)f\left(d\theta\right) & =\hat{\lambda}_{t}^{f}\cdot\mathcal{N}\left(\theta;\mu_{t}^{f},\left(\Sigma_{\mathrm{pop}}^{-1}+S_{t}^{R}\right)^{-1}\right)d\theta,\nonumber \\
\hat{\lambda}_{t}^{f} & =\int\hat{\lambda}_{t}\left(h,\theta,R\right)f\left(d\theta\right)\nonumber \\
 & =h\sqrt{\frac{\left|Z_{t}^{R}\right|}{\left|R\right|}}\exp\left(-\frac{1}{2}\left\Vert c-H\mu_{t}\right\Vert _{Z_{t}^{R}}^{2}\right)\label{eq:total-rate-gauss-gauss}
\end{align}
where
\begin{align}
Z_{t}^{R} & \triangleq\left(\Sigma_{\mathrm{pop}}+\left(S_{t}^{R}\right)^{-1}\right)^{-1}=\left(\Sigma_{\mathrm{pop}}+R^{-1}+H\Sigma_{t}H^{T}\right)^{-1}\label{eq:Z}\\
\mu_{t}^{f} & \triangleq\left(\Sigma_{\mathrm{pop}}^{-1}+S_{t}^{R}\right)^{-1}\left(S_{t}^{R}\mu_{t}+\Sigma_{\mathrm{pop}}^{-1}c\right).\nonumber
\end{align}
Substituting into (\ref{eq:mean-c-gaussian})-(\ref{eq:var-c-gaussian})
and simplifying yields the following continuous update equations
\begin{align}
d\mu_{t}^{\mathrm{c}} & =\Sigma_{t}H^{T}Z_{t}^{R}\left(H\mu_{t}-c\right)\hat{\lambda}_{t}^{f}dt\label{eq:mean-c-gauss-gauss}\\
d\Sigma_{t}^{\mathrm{c}} & =\Sigma_{t}H^{T}\Big(Z_{t}^{R}-Z_{t}^{R}\left(H\mu_{t}-c\right)\left(H\mu_{t}-c\right)^{T}Z_{t}^{R}\Big)H\Sigma_{t}\hat{\lambda}_{t}^{f}dt,\label{eq:var-c-gauss-gauss}
\end{align}
where $Z_{t}^{R}$ and $\hat{\lambda}_{t}^{f}$ are given by (\ref{eq:Z})
and (\ref{eq:total-rate-gauss-gauss}), respectively. These updates
generalize the single-sensor updates (\ref{eq:mean-c-single})-(\ref{eq:var-c-single}),
with the population center $c$ taking the place of the single location
parameter, and $Z_{t}^{R}$ substituting $S_{t}^{R}$. The single
sensor case is obtained when $\Sigma_{\mathrm{pop}}=0$.

It is illustrative to consider these equations in the scalar case
$m=n=1$, with $H=1$. Letting $\sigma_{t}^{2}=\Sigma_{t},\sigma_{\mathrm{r}}^{2}=R^{-1},\sigma_{\mathrm{pop}}^{2}=\Sigma_{\mathrm{pop}}$
yields
\begin{align}
d\mu_{t}^{\mathrm{c}} & =\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}+\sigma_{\mathrm{pop}}^{2}}\left(\mu_{t}-c\right)\hat{\lambda}_{t}^{f}dt,\label{eq:mean-c-gauss-gauss-1d}\\
d\sigma_{t}^{2,\mathrm{c}} & =\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}+\sigma_{\mathrm{pop}}^{2}}\left(1-\frac{\left(\mu_{t}-c\right)^{2}}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}+\sigma_{\mathrm{pop}}^{2}}\right)\sigma_{t}^{2}\hat{\lambda}_{t}^{f}dt,\label{eq:var-c-gauss-gauss-1d}
\end{align}
where
\[
\hat{\lambda}_{t}^{f}=h\sqrt{2\pi\sigma_{\mathrm{r}}^{2}}\mathcal{N}\left(\mu_{t};c,\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}+\sigma_{\mathrm{pop}}^{2}\right).
\]

Figure \ref{gaussian} demonstrates the continuous update terms (\ref{eq:mean-c-gauss-gauss-1d})-(\ref{eq:var-c-gauss-gauss-1d})
as a function of the current mean estimate $\mu_{t}$, for various
values of the population variance $\sigma_{\mathrm{pop}}^{2}$, including
the case of a single sensor, $\sigma_{\mathrm{pop}}^{2}=0$. The continuous
update term $d\mu_{t}^{\mathrm{c}}$ pushes the posterior mean $\mu_{t}$
away from the the population center $c$ in the absence of spikes.
This effect weakens as $\left|\mu_{t}-c\right|$ grows due to the
factor $\hat{\lambda}_{t}^{f}$, consistent with the idea that far
from $c$, the lack of events is less surprising, hence less informative.
The continuous variance update term $d\sigma_{t}^{2,\mathrm{c}}$
increases the variance when $\mu_{t}$ is near $\theta$, otherwise
decreases it. This stands in contrast with the Kalman-Bucy filter,
where the posterior variance cannot increase when estimating a static
state.

\begin{figure*}
\subfloat[Gaussian population ((\ref{eq:mean-c-gauss-gauss-1d})-(\ref{eq:var-c-gauss-gauss-1d})).
Parameters are $c=0,\sigma_{\mathrm{r}}^{2}=0.25,h=1,\sigma_{t}^{2}=1$.
The case $\sigma_{\mathrm{pop}}^{2}=0$ corresponds to a single sensor.]{\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.45\textwidth]{continuous-terms-gaussian}\label{gaussian}}\hfill{}\subfloat[{Uniform population on an interval ((\ref{eq:mean-c-interval})-(\ref{eq:var-c-interval})).
Parameters are $\left[a,b\right]=\left[-1,1\right],h=1,\sigma_{t}^{2}=0.01$.}]{\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.45\textwidth]{continuous-terms-uniform}

\label{interval}}\caption{Continuous update terms as a function of the current posterior mean
estimate, for a 1-d state observed through a population of Gaussian
sensors (\ref{eq:guassian-sensor}). The population density is Gaussian
on the left plot, and uniform on the interval $\left[-1,1\right]$
on the right plot. The bottom plots shows the population density $f\left(d\theta\right)/d\theta$
and a tuning curve $\lambda_{t}\left(x;\theta\right)$ for $\theta=0$.}
\label{continuous-terms}
\end{figure*}

\subsubsection{Uniform population on an interval}

In this case $R,h$ are constant, and $f\left(d\theta\right)$ is
given in (\ref{eq:interval}). Here we assume $m=n=1$ and $H=1$,
so (\ref{eq:mean-c-gaussian})-(\ref{eq:var-c-gaussian}) take the
form
\begin{align*}
d\mu_{t}^{\mathrm{c}} & =\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}}\int_{\theta\in\mathbb{R}}\left(\mu_{t}-\theta\right)\hat{\lambda}_{t}\left(\theta\right)f\left(d\theta\right)dt\\
d\sigma_{t}^{2,\mathrm{c}} & =\left(\hat{\lambda}_{t}^{f}-\frac{\int\left(\theta-\mu_{t}\right)^{2}\hat{\lambda}_{t}\left(\theta\right)f\left(d\theta\right)}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}}\right)\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}}\sigma_{t}^{2}dt\\
\hat{\lambda}_{t}\left(\theta\right) & =h\sqrt{2\pi\sigma_{\mathrm{r}}^{2}}\mathcal{N}\left(\theta;\mu_{t},\sigma_{\mathrm{r}}^{2}+\sigma_{t}^{2}\right),
\end{align*}
where $\sigma_{t}^{2}=\Sigma_{t},\sigma_{\mathrm{r}}^{2}=R^{-1}$,
and we suppressed the dependence of $\hat{\lambda}$ on $h,R$ from
the notation, since $h,R$ are fixed. Let $\phi\left(x\right)=\mathcal{N}\left(x;0,1\right),\Phi\left(x\right)=\int_{-\infty}^{x}\phi$.
A straightforward computation using the moments of the truncated Gaussian
distribution yields

\begin{align*}
\hat{\lambda}_{t}^{f}=\int_{a}^{b}\hat{\lambda}_{t}\left(\theta\right)d\theta & =h\sqrt{2\pi\sigma_{\mathrm{r}}^{2}}Z_{t},\\
\int\left(\theta-\mu_{t}\right)\hat{\lambda}_{t}\left(\theta\right)f\left(d\theta\right) & =-h\sqrt{2\pi\sigma_{\mathrm{r}}^{2}}z_{t}\sqrt{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}},\\
\int\left(\theta-\mu_{t}\right)^{2}\hat{\lambda}_{t}\left(\theta\right)f\left(d\theta\right) & =\left[\hat{\lambda}_{t}^{f}-h\sqrt{2\pi\sigma_{\mathrm{r}}^{2}}z'_{t}\right]\left(\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}\right),
\end{align*}
where
\begin{align*}
\alpha_{t} & =\frac{a-\mu_{t}}{\sqrt{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}}}\,, & Z_{t} & =\Phi\left(\beta_{t}\right)-\Phi\left(\alpha_{t}\right),\\
\beta_{t} & =\frac{b-\mu_{t}}{\sqrt{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}}}\,, & z_{t} & =\phi\left(\beta_{t}\right)-\phi\left(\alpha_{t}\right),\\
 &  & z'_{t} & =\beta_{t}\phi\left(\beta_{t}\right)-\alpha_{t}\phi\left(\alpha_{t}\right),
\end{align*}
yielding
\begin{align}
d\mu_{t}^{\mathrm{c}} & =h\sqrt{2\pi\sigma_{\mathrm{r}}^{2}}\sqrt{\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}}}z_{t}\sigma_{t}dt\label{eq:mean-c-interval}\\
d\sigma_{t}^{2,\mathrm{c}} & =h\sqrt{2\pi\sigma_{\mathrm{r}}^{2}}\frac{\sigma_{t}^{2}}{\sigma_{t}^{2}+\sigma_{\mathrm{r}}^{2}}z'_{t}\sigma_{t}^{2}dt\label{eq:var-c-interval}
\end{align}

Figure \ref{interval} demonstrates the continuous update terms (\ref{eq:mean-c-interval})-(\ref{eq:var-c-interval})
as a function of the current mean estimate $\mu_{t}$. When the mean
estimate is around an endpoint of the interval, the mean update $\mu_{t}^{\mathrm{c}}$
pushes the posterior mean outside the interval in the absence of spikes.
The posterior variance $\sigma_{t}^{2}$ decreases outside the interval,
where the absence of spikes is expected, and increases inside the
interval, where it is unexpected\footnote{This holds only approximately, when the tuning curve width is not
too large relative to the size of the interval. For wider tuning curves
the behavior becomes similar to the single sensor case.}. When the posterior mean is not near the interval endpoints, the
updates are near zero, consistently with the uniform population case
(\ref{eq:mean-var-uniform}).

\subsubsection{Finite Mixtures}

Assume $f$ is any finite mixture of measures
\[
f\left(dy\right)=\sum_{i}\alpha_{i}f_{i}\left(dy\right),
\]
where $f_{i}$ is of one of the above forms. We note that the continuous
updates (\ref{eq:mean-c})-(\ref{eq:var-c}) are linear in $f$, so
they are obtained by the appropriate weighted sums of the filters
derived above for the various special forms of $f_{i}$.

In particular, for a general mixture of the first three forms\footnote{In the one-dimensional case, the fourth form may be included similarly},
\begin{align*}
f\left(dh,dR,d\theta\right) & =\sum_{i}\alpha_{i}^{\delta}\delta_{h_{i}}\left(dh\right)\delta_{R_{i}}\left(dR\right)\delta_{\theta_{i}}\left(d\theta\right)\\
 & \quad+\sum_{i}\alpha_{i}^{U}\delta_{h_{i}^{U}}\left(dh\right)\delta_{R_{i}^{U}}\left(dR\right)d\theta\\
 & \quad+\sum_{i}\alpha_{i}^{\mc N}\delta_{h_{i}^{\mc N}}\left(dh\right)\delta_{R_{i}^{\mc N}}\left(dR\right)\mathcal{N}\left(\theta;c_{i},\Sigma_{\mathrm{pop}}^{\left(i\right)}\right)d\theta,
\end{align*}
the resulting continuous update terms are given by
\begin{align*}
d\mu_{t}^{\mathrm{c}}=\Sigma_{t}H^{T}\Bigg( & \sum_{i}\alpha_{i}^{\delta}S_{t}^{i}\left(H\mu_{t}-\theta_{i}\right)\hat{\lambda}_{t}\left(h_{i},\theta_{i},R_{i}\right)\\
\qquad+ & \sum_{i}\alpha_{i}^{\mc N}Z_{t}^{i}\left(H\mu_{t}-c_{i}\right)\hat{\lambda}_{t}^{f,i}\bigg)dt,
\end{align*}
\begin{align*}
d\Sigma_{t}^{\mathrm{c}}=\Sigma_{t}H^{T}\Bigg\{ & \sum_{i}\alpha_{i}^{\delta}\left[S_{t}^{i}-S_{t}^{i}\left(H\mu_{t}-\theta_{i}\right)\left(H\mu_{t}-\theta_{i}\right)^{T}S_{t}^{i}\right]\hat{\lambda}_{t}\left(h_{i},\theta_{i},R_{i}\right)\\
 & +\sum_{i}\alpha_{i}^{\mc N}\left[Z_{t}^{i}-Z_{t}^{i}\left(H\mu_{t}-c_{i}\right)\left(H\mu_{t}-c_{i}\right)^{T}Z_{t}^{i}\right]\hat{\lambda}_{t}^{f,i}\Bigg\} H\Sigma_{t}dt,
\end{align*}
where $\hat{\lambda}_{t}\left(h,\theta,R\right)$ is given by (\ref{eq:self-rate-gaussian}),
$S_{t}^{i}=S_{t}^{R_{i}},Z_{t}^{i}=Z_{t}^{R_{i}^{\mc N}}$ as defined
in (\ref{eq:S}),(\ref{eq:Z}), and similarly,
\[
\hat{\lambda}_{t}^{f,i}=h_{i}^{\mc N}\sqrt{\frac{\left|Z_{t}^{i}\right|}{\left|R_{i}^{\mc N}\right|}}\exp\left(-\frac{1}{2}\left\Vert c_{i}-H\mu_{t}\right\Vert _{Z_{t}^{i}}^{2}\right)
\]
(cf. (\ref{eq:total-rate-gauss-gauss})). The uniform components of
the mixture do not appear explicitly in the filter, since the continuous
update part vanishes for the uniform case (see (\ref{eq:mean-var-uniform})),
though they would influence the discontinuous update terms through
their effect on the statistics of $N$.

\section{Evaluation of filter}

\subsection{Examples and comparison to particle filter}

\begin{figure*}
\subfloat[high rate of events: $h=1000$]{\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.45\textwidth]{pf_rate-1000_4,0}\label{high-rate}}\hfill{}\subfloat[low rate of events: $h=2$]{\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.45\textwidth]{pf_rate-2_0,0}\label{low-rate}}\caption{Two examples of a linear one-dimensional process observed through
a Gaussian population (\ref{eq:gaussian-f}) of Gaussian sensors (\ref{eq:guassian-sensor}),
filtered using the ADF approximation (namely, equations (\ref{eq:mean-p})-(\ref{eq:var-p}),(\ref{eq:mean-c-gaussian})-(\ref{eq:var-c-gaussian}),(\ref{eq:mean-c-gauss-gauss})-(\ref{eq:var-c-gauss-gauss})),
and using a particle filter. Each dot correspond to a spike with the
vertical location indicating the sensor's location parameter $\theta$.
The approximate posterior means obtained from ADF and particle filtering
are shown in blue and red respectively, with the corresponding posterior
standard deviations in shaded areas of the respective colors. The
curves to the right of the graph show the preferred stimulus density
(black), and a sensor's response rate function (\ref{eq:guassian-sensor})
centered at $\theta=0$ (gray), arbitrarily normalized to the same
height for visualization. The bottom graph shows the posterior variance.
Parameters used in both examples: $a=-0.1,d=1,H=1,\sigma_{\mathrm{pop}}^{2}=4,R^{-1}=0.25,c=0,\mu_{0}=0,\sigma_{0}^{2}=1$
(note the different extent of the time axis). The observed processes
were initialized from their steady-state distribution. The dynamics
were discretized with time step $\Delta t=10^{-3}$. The particle
filter uses 1000 particles with systematic resampling (see, e.g.,
\cite{DouJoh09}) at each time step.}

\label{1d-filtering}
\end{figure*}

Since the filter (\ref{eq:mean-p})-(\ref{eq:var-N}) is based on
an assumed density approximation, its results may be inexact. We tested
the accuracy of the filter in the Gaussian population case (\ref{eq:mean-c-gauss-gauss})-(\ref{eq:var-c-gauss-gauss}),
by numerical comparison with Particle Filtering (PF) \cite{DouJoh09}.

Figure \ref{1d-filtering} shows two examples of filtering a one-dimensional
process observed through a Gaussian population of Gaussian sensors
(\ref{eq:gaussian-f}), using both the ADF approximation (\ref{eq:mean-c-gauss-gauss-1d})-(\ref{eq:var-c-gauss-gauss-1d})
and a Particle Filter (PF) for comparison. See the figure caption
for precise details.

Figure \ref{ks} shows the distribution of approximation errors and
their relation to the deviation of the posterior from Gaussian. The
deviation of the particle distribution from Gaussian is quantified
using the Kolmogorov-Smirnov (KS) statistic $\sup_{x}\left|F\left(x\right)-G\left(x\right)\right|$
where $F$ is the particle distribution cdf and $G$ is the cdf of
a Gaussian matching $F$'s first two moments. The approximation errors
plotted are the relative error in the mean estimate $\epsilon_{\mu}\triangleq\left(\mu_{\mathrm{ADF}}-\mu_{\mathrm{PF}}\right)/\sigma_{\mathrm{PF}}$,
and the error in the posterior standard deviation estimate $\epsilon_{\sigma}\triangleq\sigma_{\mathrm{ADF}}/\sigma_{\mathrm{PF}}$,
where $\mu_{\mathrm{ADF}},\mu_{\mathrm{PF}},\sigma_{\mathrm{ADF}},\sigma_{\mathrm{PF}}$
are, respectively, the posterior mean obtained from ADF and PF, and
the posterior variance obtained from ADF and PF. The observed mean
and standard deviation of the estimation errors are $0.0018\pm0.0989$
for $\epsilon_{\mu}$ and $1.010\pm0.101$ for $\epsilon_{\sigma}$.
The results suggest that the largest errors occur in rare cases where
the posterior diverges significantly from Gaussian, and involve an
overestimation of the posterior variance. Similar results were obtained
with different parameters.

\begin{figure}
\includegraphics[width=1\columnwidth]{KS-errs-1000_0}\caption{Approximation errors (relative to particle filter) vs. KS statistic
of particle distribution. The KS statistic is plotted against the
estimation errors $\epsilon_{\mu},\epsilon_{\sigma}$ (see main text).
Each point corresponds to a single simulation time step. The histograms
show the distribution of the KS statistic, of $\epsilon_{\mu}$ and
of $\epsilon_{\sigma}$. Results obtained from 100 trials of length
$T=1$, with parameters as in Figure \ref{high-rate}. }
\label{ks}
\end{figure}

\subsection{Information gained between spikes}

The filters derived above for various population distributions differ
only in the continuous update terms, which modify the posterior between
spikes beyond the prior terms derived from the state dynamics. We
may therefore interpret this term as corresponding to information
gained by the absence of spikes. The continuous update term vanishes
in the uniform population filter of \cite{RhoSny1977} (see (\ref{eq:mean-var-uniform})),
so that in the uniform case, between spikes the posterior dynamics
are identical to the prior dynamics. This reflects the fact that lack
of spikes in a time interval is an indication that the total firing
rate is low; in the uniform population case, this is not informative,
since the total firing rate is independent of the state.

Figure \ref{adf-vs-uniform} (left) illustrates the contribution of
the continuous update terms (\ref{eq:mean-c-gauss-gauss})-(\ref{eq:var-c-gauss-gauss})
to filter performance. A static scalar state is observed by a Gaussian
population (\ref{eq:guassian-sensor})-(\ref{eq:gaussian-f}), and
filtered twice: once with the correct value of $\sigma_{\mathrm{pop}}^{2}=\Sigma_{\mathrm{pop}}$,
and once with $\sigma_{\mathrm{pop}}\rightarrow\infty$, which yields
$\mu_{t}^{\mathrm{c}}=0$, recovering the uniform population filter
of \cite{RhoSny1977}. Between spikes, the ADF estimate moves away
from the population center $c=0$, whereas the uniform coding estimate
remains fixed. The size of this effect decreases with time, as the
posterior variance estimate (not shown) decreases. The reduction in
filtering errors gained from the continuous update terms is illustrated
in Figure \ref{adf-vs-uniform} (right). Despite the approximation
involved, the full filter significantly outperforms the uniform population
filter. The difference disappears as $\sigma_{\mathrm{pop}}$ increases
and the population becomes uniform.

\begin{figure*}
\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.4\textwidth]{uniform-vs-adf}\hfill{}\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.4\textwidth]{ADF-uniform-MSE-1000}\caption{\textbf{Left} Illustration of information gained between spikes. A
static state $X_{t}=0.5$, shown in a dotted line, is observed by
a Gaussian population (\ref{eq:gaussian-f}) of Gaussian sensors (\ref{eq:guassian-sensor}),
and filtered twice: with the correct value $\sigma_{\mathrm{pop}}^{2}=0.5$
(``ADF'', solid blue line), and with $\sigma_{\mathrm{pop}}^{2}=\infty$
(``Uniform coding filter'', dashed line), which is equivalent to
ignoring the continuous update terms. Both filters are applied to
the same random realization of the observation process. The curves
to the right of each graph show the preferred stimulus density (black),
and a tuning curve centered at $\theta=0$ (gray). Both filters are
initialized with $\mu_{0}=0,\sigma_{0}^{2}=1$. \textbf{Right} Comparison
of MSE for the ADF filter and the uniform coding filter. The vertical
axis shows the integral of the square error integrated over the time
interval $\left[5,10\right]$, averaged over 1000 trials. Shaded areas
indicate estimated errors, computed as the sample standard deviation
divided by the square root of the number of trials. Parameters in
both plots are $a=d=0,c=0,\sigma_{\mathrm{pop}}^{2}=0.5,\sigma_{\mathrm{r}}^{2}=0.1,H=1,h=10$.}
\label{adf-vs-uniform}
\end{figure*}

\section{Encoding}

\subsection{Motivation}

We demonstrate the use of the Assumed Density Filter in determining
optimal encoding strategies, i.e., selecting the optimal sensor population
parameters $\phi$ (see Section \ref{sub:Encoding-and-decoding} above).
The study of optimal encoding has both biological and engineering
motivations. In the context of neuroscience, the optimal encoding
strategy may serve as a model for the observed characteristics of
sensory neurons, or for observed sensory adaptation (e.g., \cite{DayAbb05}).
Several works have studied optimal neural encoding using Fisher information
as a performance criterion. Fisher information is easily computed
from tuning functions in the case of static state (e.g., \cite{DayAbb05}),
and this approach completely circumvents the need to model the decoding
process. For example, in \cite{GanSim14}, the authors consider a
scalar static state observed by a parameterized population of neurons,
and analytically optimize the parameters based on Fisher information
under a constraint on the total expected rate of spikes. A few works
(e.g. \cite{YaeMei10,Susemihl2014}) consider the more difficult problem
of direct minimization of the MSE rather than the Fisher information,
by solving the corresponding filtering problem and measuring the MSE
of different encoding parameters using Monte Carlo (MC) simulations.
The filtering problem is made tractable in \cite{YaeMei10} and \cite{Susemihl2014}
by assuming a uniform population (equivalent to case 2 in Section
\ref{sub:special-cases}) above. However, sensory populations are
often non-uniform (e.g. \cite{Brand2002}), and sensory adaptation
often modifies tuning curves non-uniformly (e.g. \cite{Benucci2013}),
which motivates our more general setting. Optimization of sensor configuration
and coding is also increasingly studied in engineering contexts (e.g.,
\cite{Yuksel2013,andrievsky2010,Mourikis2006,Ny2011}), as networked
filtering and control are becoming ubiquitous. The latter studies
are usually concerned with continuous observations (often linear),
rather than PP based observations using heterogeneous biologically
motivated tuning functions as is done here.

\subsection{Encoding example}

To illustrate the use of ADF for the encoding problem, we consider
a simple example using a Gaussian population (\ref{eq:gaussian-f}).
We will study optimal encoding issues in more detail in a sequel paper.

Previous work using a finite neuron population and a Fisher information-based
criterion \cite{HarMcAlp04} has suggested that the optimal distribution
of preferred stimuli depends on the prior variance. When it is small
relative to the tuning curve width, optimal encoding is achieved by
placing all preferred stimuli at a fixed distance from the prior mean.
On the other hand, when the prior variance is large relative to the
tuning curve width, optimal encoding is uniform (see figure 2 in \cite{HarMcAlp04}).
These results are consistent with biological observations reported
in \cite{Brand2002} concerning the encoding of aural stimuli.

Similar results are obtained with our model, as shown in Figure \ref{harper}.
Whereas \cite{HarMcAlp04} implicitly assumed a static state in the
computation of Fisher information, we use a time-varying scalar state.
The state obeys the dynamics
\[
dX_{t}=aX_{t}dt+d\,dW_{t}\quad\left(a<0\right),
\]
and is observed through a Gaussian population (\ref{eq:gaussian-f})
and filtered using the ADF approximation. In this case, optimal encoding
is interpreted as the simultaneous optimization of the population
center $c$ and the population variance $\Sigma_{\mathrm{pop}}$.
The process is initialized so that it has a constant prior distribution,
its variance given by $d^{2}/\left(2\left|a\right|\right)$. In Figure
\ref{harper} (left), the steady-state prior distribution is narrow
relative to the tuning curve width, leading to an optimal population
with a narrow population distribution far from the origin. In Figure
\ref{harper} (right), the prior is wide relative to the tuning curve
width, leading to an optimal population with variance that roughly
matches the prior variance.

Our approach, though more computationally expensive, offers two advantages
over the Fisher information-based method which is used in \cite{HarMcAlp04}
and which is prevalent in computational neuroscience. First, the simple
computation of Fisher information from tuning curves, commonly used
in the neuroscience literature, is based on the assumption of a static
state, whereas our method can be applied in a fully dynamic context,
including the presence of observation-dependent feedback. Second,
our approach allows the minimization of arbitrary criteria, including
the direct minimization of posterior variance or Mean Square Error
(MSE). Although, under appropriate conditions, Fisher information
approaches the MSE in the limit of infinite decoding time, it may
be a poor proxy for the MSE for finite decoding times (e.g., \cite{Bethge2002,YaeMei10}),
which are of particular importance in natural settings and in control
problems.

\begin{figure*}
\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.55\textwidth]{Posterior_variance,_narrow_prior}\hfill{}\includegraphics[bb=0bp 0bp 534bp 390bp,width=0.55\textwidth]{Posterior_variance,_wide_prior}\caption{Optimal population distribution depends on prior variance relative
to tuning curve width. A scalar state with dynamics $dX_{t}=aX_{t}+d\,dW_{t}$
($a=-0.05$) is filtered with tuning curves parameters $h=50,\sigma_{\mathrm{r}}=1$
and preferred stimulus density $\mathcal{N}(c,\sigma_{\mathrm{pop}}^{2})$.
The process is initialized from its steady state distribution, $\protect\mc N\left(0,d^{2}/\left(-2a\right)\right)$.
Both graphs show the posterior standard deviation derived from the
ADF approximation, relative to the prior standard deviation $\sigma_{\mathrm{p}}$.
In the left graph, $d=0.1$ so that the prior variance is $0.1$,
whereas on the right, $d=1$, so that the prior variance is $10$.
In both cases the filter is initialized with the correct prior, and
the posterior variance is averaged over the time interval $\left[1,2\right]$
and across 1000 trials for each data point. Only non-negative values
of $c$ were simulated, but note that the process is symmetric about
zero, so that the full plots would also be symmetric about $c=0$.
The areas colored in white in the right plot correspond to parameters
where the computed posterior variance exceeded the prior variance.
This is due to poor performance of the ADF approximation for these
parameter values, in cases where no spikes occur and the true posterior
becomes bimodal.}
\label{harper}
\end{figure*}

\section{Conclusions}

We have introduced an analytically tractable approximation to point
process filtering, allowing us to gain insight into the generally
intractable infinite-dimensional filtering problem. The approach enables
the derivation of near-optimal encoding schemes going beyond the previously
studied case of uniform population. The framework is presented in
continuous time, circumventing temporal discretization errors and
numerical imprecision in sampling-based methods, applies to fully
dynamic setups, and directly estimates the MSE rather than lower bounds
to it. It successfully explains observed experimental results, and
opens the door to many future predictions. Moreover, the proposed
strategy may lead to practically useful decoding of spike trains.

\bibliographystyle{IEEEtran}

redpajama_set_name: RedPajamaArXiv


text:
Voltaren makers to pay $4.5m for claiming Osteo Gel better than Emulgel
This was published 7 months ago
By Matt Bungard
Updated May 28, 2020 — 6.58pm first published at 6.02pm
The pharmaceutical giants who produce Voltaren sports gel have been ordered by a Federal Court to pay $4.5 million in damages after they were ruled to have breached Australian Consumer Law.
GlaxoSmithKlein and Novartis were taken to court over claims that Voltaren's Osteo Gel was more effective in treating pain than Voltaren Emulgel, while costing 33 per cent more.
Voltaren Osteo Gel and Emulgel were found to contain the same ingredients.
From January 2012 to March 2017, their Osteo Gel product was marketed as being specifically formulated and more effective than the cheaper Emulgel in treating osteoarthritis-related pain and inflammation, even though both had the same active ingredients.
"It's really disappointing that post-Nurofen, there is still behaviour where we have products with the same active ingredients stating two different abilities and two different remedies," ACCC chairman Rod Sims said at the time.
In 2016, pharmaceutical company Reckitt Benckiser, which manufactures Nurofen, was ordered to pay a $6 million penalty for misleading consumers with its specific pain relief range, which was confirmed after a final appeal to the High Court was rejected in August of 2017.
The $6 million corporate penalty remains the highest-ever contested fine under Australian Consumer Law.
The highest penalty the Federal Court has ordered for breaches of Australian Consumer Law were $10 million, in two separate ACCC actions against Coles in 2014, and Ford in 2018. They have since been passed by Volkswagen ($125 million), Empower Institute ($26.5 million) and We Buy Houses and its sole director Rick Otton ($18 million).
The highest penalty the Federal Court has ordered for breaches of competition law is $46 million in an ACCC action against Yazaki.
Both Osteo Gel and Emulgel products contain the same active ingredient: diclofenac dimethylammonium gel 11.6mg/g, which acts in a non-specific manner to reduce local pain and inflammation wherever it is applied.
"Novartis and GSK's claims, which made Osteo Gel seem more effective than Emulgel, were unacceptable," ACCC Commissioner Sarah Court said on Thursday.
She said the claims were particularly concerning because consumers were potentially misled into paying more for an identical product believing it was more effective.
"Voltaren Osteo Gel and Voltaren Emulgel were essentially the same gel and were equally effective in treating osteoarthritis symptoms."
The misleading claims were made on product packaging, as well as on the Voltaren website by both companies.
Novartis also made the claims on the My Joint Health website.
The Arthritis Australia website describes osteoarthritis as a condition that affects the whole joint including bone, cartilage, ligaments and muscles and although it is often described as "wear and tear", the condition is now thought to be the result of a joint working extra hard to repair itself.
"We believe this case highlights the importance of ensuring that marketing claims for therapeutic products are legitimate and evidence-based, and that vulnerable consumers are not misled into paying a price premium for a product which offers no additional benefit," an Arthritis Australia spokesperson told the Herald.
Matt Bungard
Matt Bungard is a journalist at The Sydney Morning Herald.

redpajama_set_name: RedPajamaCommonCrawl


text:
News Channel 3 TV Schedule Autos
State grants development money to Tennessee cities, counties
Posted 5:40 pm, October 11, 2019, by Peter Fleischer, Updated at 06:38PM, October 14, 2019
*** This story is an updated version with multiple corrections after our first version of the story contained various inaccuracies. We regret the errors, as accuracy is of the utmost importance to WREG News Channel 3 ***
MEMPHIS, Tenn. — Gov. Bill Lee approved millions of dollars in grant money Friday morning that will help counties around the state with infrastructure improvements and safety initiatives. But when the final list came out, much of West Tennessee was left off.
More than $26 million was given out to 66 different counties in the state.
Community leaders agreed, they'll take all the help they can get, but there are different routes for different types of areas.
More than two-thirds of the counties across Tennessee received a slice of the almost $30 million that Gov. Lee approved Friday morning. In West Tennessee, Lauderdale County, Brighton, Gallaway, Stanton, Bolivar, Middleton, Toone, Adamsville and Michie all received grant money.
Some grants did go to cities and towns in the area, including about half a million dollars for a sewer project in Bolivar.
In total, West Tennessee will receive about $8.5 million, which does not include Memphis or Shelby County.
WREG reached out to the Tennessee Department of Economic and Community Development to understand why so little of that money made its way to the Mid-South. An official explained over the phone there are "entitlement cities" that are not eligible for the community development block grants because they already receive funding from the Housing and Urban Development Department.
That means Memphis and Shelby County aren't considered for the state grants. State Rep. Antonio Parkinson explained the difficulties of balancing funding, saying avoiding the state process can actually be beneficial.
"The good thing about it is, we don't have to go through the same process that some of the smaller counties have to go through," Rep. Parkinson said. "So we get our monies directly from the federal government without having the state government as a middle man."
Through the federal program, Memphis is receiving $6.6 million, and Shelby County is receiving $1.1 million.
A state official told WREG they balance the decision by combining community need and project impact, terms that at least this year have left a large portion of southwest corner of the state out to dry. If anything, it provides a clear window in the competitive state process, where everyone wants a piece of the pie.
"If you're fighting for something that's going to shore up a weak spot, that weak spot could be something that's pulling in time and budget from the state now," Rep. Parkinson said.
Topics: development grant, Gov. Bill Lee
News, weather, and sports headlines from our team
Two children killed in Hickory Hill shooting identified as Collierville Schools students
'Jeopardy!' and 'Inside Edition' to air later tonight
Gov. Lee decides Tennessee will continue accepting refugees
Tennessee governor says state will keep resettling refugees
Gov. Lee wants to keep 'best and brightest' teachers in Tennessee classrooms
Lee remains hopeful on school voucher timeline
Tennessee governor working to dip into welfare reserve
First of many federal broadband grants announced in Tennessee
Tennessee gov appoints former GOP lawmaker to parole board
Tennessee Medicaid block grant plan heads for federal review
Think tank finds $732 million for Tennessee's poor is unused
Tennessee governor says block grant opponents 'misinformed'
Lee announces rural broadband expansion in West Tennessee
Tennessee doling out $5.9M in grants to 117 charter schools
Tennessee congressman seeks answers on welfare surplus


Anyone for Heaven? By John Blanchard and published by Evangelical Press.
As a friend of mine often reminds people, 'We are not here to stay; we are here to go.' I am writing these words in a Consultant's waiting room in a London hospital. Some of the other patients are showing signs of great physical distress, moving reminders that the whole world is a hospital in which everyone is a terminal patient.
The South African surgeon Christiaan Barnard performed the world's first human heart transplant operation in 1967 and the first double-heart transplant in 1974. Later in life he said, 'At 79, people ask me, "Where do you go from here?" I say to them, "I'm on the waiting list. I don't exactly know where I am on that, or where I'm going. But I'm on it".' So are you! Even as you read these words you are fighting a losing battle with death, whatever your age or present state of health. What will happen to you when death wins..?


text:
How the Aquila brand changed the way hotels buy hotels
When a brand can get away with making such a drastic change in the way people shop, it's hard to believe that the trend is so widespread.
In fact, it was only in the last few years that hotels started using the Aquilas brand, and it's only now that brands like Marriott and Hilton are catching on.
Now, hotels are taking the Aquiles brand to the next level, using them to make a wide variety of items from toiletries to water purifiers to bathtubs.
Here are 10 of the coolest hotels using Aquila to sell their products directly.1.
Aquila hotel supplies The Aquila Hotel Supply Co. has made a name for itself selling direct-to-consumer hotels in the U.S. through its website.
Its products include water purification systems and toilet paper, as well as laundry detergent and detergent wipes.
Aquilases soap and shampoo, and the brand's line of bath products are also available for direct purchase on its website, though customers will have to purchase the product at a local retailer, and can't get them directly from Aquila.
The brand's "Made in the USA" tag makes it easy to identify products as being made in the United States.2.
Aquils bath towels AquilAs brands, hotels tend to follow trends, but the Aquils line has taken a more radical approach to selling its products online.
Instead of simply copying the Aquilla brand's name, Aquilaks products have been designed with the help of designers from other brands and companies.
This is the first time Aquilays brand has been used to sell directly, and its name alone is enough to make some customers say, "wow."3.
Aqualife Aquilife has a reputation for being a great brand, but its website makes it clear that it is not the only brand in the market.
Aquilon is a company that sells direct-from-the-brand products, but they don't look like the usual Aquila products.
Instead, Aquilon uses a "Made by" tag, with the brand name on the front.
And while Aquilay's brands are usually branded by brands that are owned by major corporations, the Aquilon name stands out.
"Made for you, Aquilan, the future of Aquila," the company's website reads.4.
Aquileam Aquilaseram is a boutique hotel chain that has expanded into several other parts of the world.
Its website offers many of its hotels as direct-sell items, though they also sell products that are part of Aquileas brand.
Some of the products, like the AquaLite, have the brand on the back.
But it also offers products made by other companies that Aquila is affiliated with, such as the AquaServe line and the Aqualiks.5.
Aquilia hotels and resorts hotels in a hotel room, but in the real world, the hotel room is a much bigger part of the experience.
In the real-world, hotels and motels are usually a part of homes, so if a hotel or resort has a lot of guests, it may be a good idea to sell a lot more of the hotel or motel product than just the items that you buy directly from the hotel.
So, some hotels may sell a few products directly to customers.
In that case, some of the items may be "supplied directly by the hotel," which is an industry term for an item that is directly manufactured by the company.
Some hotels have a different system of sourcing products, and may only buy items from the suppliers they deal with directly.
Some resorts may use distributors that provide products directly from other suppliers, but other resorts may buy from other hotels or motels.6.
Aquilleas laundry detergents The Aquilades brand has made it possible for consumers to buy laundry detergenators directly from a brand that doesn't have to go through a retailer.
The products are not just used in hotels, either.
They are also sold through Aquila stores.
Aquiles soap and detergent, for example, are sold directly through Aquilashes stores.
And the brand has also expanded its range of other products, including soap, toothpaste, and shower heads.
The brands' packaging is also much more straightforward than a lot the other brands.
There are no expensive labels that are on the labels.
The packaging is printed on the side of the package, so you can easily see the brands name, logo, and logo design.7.
Aquilias spa products There are a few hotels in Japan that specialize in Spa services, but hotels in China and other parts are starting to take the Aquiliases brand seriously as well.
In those areas, hotels have been adding a few spa services, including a few that sell Aquila shampoo.
The Aquiliades brand offers soap and water that can be used to cleanse the scalp, as opposed to the
Tagged: aqua hotel supply, elemis hotel supplies, hotel supplies direct
How to get ready for a hurricane in Florida
When the sharjah airport closed, hotel prices went up


So want to see the worlds best wakeskater and wakeboarder 2017 edit? Of course you do. Daniel Grant is a absolute beast in freeriding and of course contest…he basically won everything.
Liquid Force is doing it big this year with their LTD line of wakeboard and bindings for 2015. Once again the product is only available at your local wakeboard retailers, but don't worry we should hopefully have some before too long. The whole line looks really good! What set up would you rock?
Peep the Remedy LTD, Peak LTD, Classic LTD, Tao LTD and their respected boot partners the Harley LTD, Shane LTD, Watson LTD, and Tao LTD.
There is also a pretty slick video that Liquid Force did – Check it out here!
Last updated by Rhett at November 29, 2017 .


Hah! I did have a color of the year in 2017. I didn't act much upon it at the time but I did buy some yarn: bright orange and hot pink. I think it's a great combination and one day it will be come a shawl. It feels good to know that the tradition isn't broken, as it has now become a thing.


text:
Bibliographic details for Lhamthanc
Page name: Lhamthanc
Date of last revision: 23 December 2021 16:38 UTC
Permanent URL: https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501
Citation styles for Lhamthanc
Lhamthanc. (2021, December 23). Tolkien Gateway, . Retrieved 11:40, January 29, 2023 from https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501.
"Lhamthanc." Tolkien Gateway, . 23 Dec 2021, 16:38 UTC. 29 Jan 2023, 11:40 <https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501>.
Tolkien Gateway contributors, 'Lhamthanc', Tolkien Gateway, , 23 December 2021, 16:38 UTC, <https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501> [accessed 29 January 2023]
Tolkien Gateway contributors, "Lhamthanc," Tolkien Gateway, , https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501 (accessed January 29, 2023).
Tolkien Gateway contributors. Lhamthanc [Internet]. Tolkien Gateway, ; 2021 Dec 23, 16:38 UTC [cited 2023 Jan 29]. Available from: https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501.
Lhamthanc, https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501 (last visited January 29, 2023).
title = "Lhamthanc --- Tolkien Gateway{,} ",
url = "https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501",
url = "\url{https://tolkiengateway.net/w/index.php?title=Lhamthanc&oldid=339501}",


While every effort is made to ensure that the information contained in this site is correct, Giles Finance makes no warranty express or implied as to the nature or accuracy of any such material to the extent permitted by applicable law.
Giles Finance & Consultancy Services may make changes to any content and information contained on this site, at any time without notice.
'Cookies' are small pieces of information that are stored in your browser on your computer's hard drive. Our cookies are used for the management and administration of the site only. They do not contain any personally identifying information.
Payable to the Buy to Let mortgage lender for the administration of the case. This can be added to the overall Buy to Let loan or paid upon submission of the application.
A charge for securing a certain Buy to Let mortgage rate, usually only found on Buy to Let fixed rates.This must be paid upon submission of the application.
A small percentage of Buy to Let mortgage lenders will also apply the following fee:Please see application fee scale below.
We will not charge you broker fees for arranging a Buy to Let Mortgage however broker fees will apply to all Commercial Mortgages. Some of the following lender fees may apply.These will be fully outlined in a key fact illustration issued.
We shall charge you application and broker fees for arranging a Commercial Mortgage. Because Commercial Mortgage Finance is mainly on bespoke basis, we do not have a set standard fee scale. However we promise to minimise your cost as best as we can.
Please fill out an enquiry form and we will call you back to discuss.
The legal ownership and any copyright or other intellectual property rights in the site and its design, text, graphics, software and all other material on the site belong to Giles Finance & Consultancy Services. You may only download content or view and/or print hard copies for private or non commercial use and no part of the site or its material shall otherwise be produced, transmitted or made available on a network without the prior written permission of Giles Finance & Consultancy Services.
Any framing of or linking to this site is prohibited without the prior written permission of Giles Finance & Consultancy Services.
No. Giles Finance do not sell, trade or rent your information to others.
Giles Finance & Consultancy Services is providing this web site on an 'as is' basis and makes no representations or warranties of any kind with respect to this website or its contents and disclaims all such representations and warranties.
The Buy to Let mortgage lender will use the same solicitor as the customer (the practice must have 2 or more partners).All legal fees are payable by the customer, as and when required by their solicitor.
OFFER AVALIABLE TO ENQUIRIES MADE VIA WEBSITE ONLY.
Giles Finance & Consultancy Services are committed to ensuring that your privacy is protected. This policy explains how Giles Finance & Consultancy Services use any information collected about you.
Giles Finance will charge you for aranging secured loans.
Please note that Giles Finance & Consultancy Services reserves the right to make changes to these prices without prior notice.Please fill out an enquiry form and we will call you back to discuss.
Payable to Giles Finance & Consultancy Services for arranging the actual property valuation which is required by the Buy to Let mortgage lender.Please note that a copy of this valuation is normally only sent to the Buy to Let mortgage lender. This fee is non-refundable must be paid upon submission of the application together with the application fee.
Please see application fee scale below.
We will provide you information on various buy to let products to help you arrive at a meaningful decision that suits your financial needs.
Giles Finance & Consultancy Services.
We are not a lender and do not act on behalf of any lender in providing mortgage information. The lender's terms will apply to any application you make through us.FSA does not regulate Buy to Let Mortgages,Commercial Mortgages, National Savings Accounts & Taxation advice.
There will be a maximum brokerage fee of 1% of loan amount on Buy to Let and Investment mortgage services. There will be a charge for Commercial Mortgage arrangements.
When you make an enquiry, send a quotation or join our mailing list you provide your name and contact details. This allows Giles Finance & Consultancy Services to fulfil the request and make a follow up telephone call or email. The information also allows Giles Finance & Consultancy Services to answer your enquiry or keep you up to date with market or product developments.