[pointprocesses] improved count statistics
This commit is contained in:
parent
849f1d8f0d
commit
c5d5d94d1b
49
pointprocesses/lecture/countexamples.py
Normal file
49
pointprocesses/lecture/countexamples.py
Normal file
@ -0,0 +1,49 @@
|
||||
import numpy as np
|
||||
from scipy.stats import poisson
|
||||
import matplotlib.pyplot as plt
|
||||
from plotstyle import *
|
||||
|
||||
|
||||
rate = 20.0
|
||||
trials = 10
|
||||
duration = 200.0
|
||||
|
||||
|
||||
def hompoisson(rate, trials, duration) :
|
||||
spikes = []
|
||||
for k in range(trials) :
|
||||
times = []
|
||||
t = 0.0
|
||||
while t < duration :
|
||||
t += np.random.exponential(1/rate)
|
||||
times.append( t )
|
||||
spikes.append( times )
|
||||
return spikes
|
||||
|
||||
|
||||
def plot_count_hist(ax, spikes, win, pmax):
|
||||
counts = []
|
||||
for times in spikes:
|
||||
c, _ = np.histogram(times, np.arange(0.0, duration, win))
|
||||
counts.extend(c)
|
||||
cb = np.arange(0.0, 15.5, 1.0)
|
||||
h, b = np.histogram(counts, cb, density=True)
|
||||
ax.bar(b[:-1], h, bar_fac, **fsA)
|
||||
ax.plot(cb, poisson.pmf(cb, rate*win), **lsBm)
|
||||
ax.plot(cb, poisson.pmf(cb, rate*win), **psBm)
|
||||
ax.text(0.9, 0.9, 'T=%.0fms' % (1000.0*win), transform=ax.transAxes, ha='right')
|
||||
ax.set_xlim(-0.5, 10.5)
|
||||
ax.set_ylim(0.0, pmax)
|
||||
ax.set_xticks(np.arange(0.0, 11.0, 5.0))
|
||||
ax.set_xlabel('Counts k')
|
||||
ax.set_ylabel('P(k)')
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
spikes = hompoisson(rate, trials, duration)
|
||||
fig, (ax1, ax2) = plt.subplots(1, 2)
|
||||
fig.subplots_adjust(**adjust_fs(fig, top=0.5, right=1.5))
|
||||
plot_count_hist(ax1, spikes, 0.02, 0.7)
|
||||
plot_count_hist(ax2, spikes, 0.2, 0.25)
|
||||
plt.savefig('countexamples.pdf')
|
||||
plt.close()
|
||||
107
pointprocesses/lecture/fanoexamples.py
Normal file
107
pointprocesses/lecture/fanoexamples.py
Normal file
@ -0,0 +1,107 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.ticker as mpt
|
||||
from plotstyle import *
|
||||
|
||||
|
||||
rate = 20.0
|
||||
trials = 20
|
||||
duration = 100.0
|
||||
dt = 0.001
|
||||
drate = 50.0
|
||||
tau = 0.1;
|
||||
|
||||
|
||||
def hompoisson(rate, trials, duration) :
|
||||
spikes = []
|
||||
for k in range(trials) :
|
||||
times = []
|
||||
t = 0.0
|
||||
while t < duration :
|
||||
t += np.random.exponential(1/rate)
|
||||
times.append( t )
|
||||
spikes.append( times )
|
||||
return spikes
|
||||
|
||||
|
||||
def inhompoisson(rate, trials, dt) :
|
||||
spikes = []
|
||||
p = rate*dt
|
||||
for k in range(trials) :
|
||||
x = np.random.rand(len(rate))
|
||||
times = dt*np.nonzero(x<p)[0]
|
||||
spikes.append( times )
|
||||
return spikes
|
||||
|
||||
|
||||
def pifspikes(input, trials, dt, D=0.1) :
|
||||
vreset = 0.0
|
||||
vthresh = 1.0
|
||||
tau = 1.0
|
||||
spikes = []
|
||||
for k in range(trials) :
|
||||
times = []
|
||||
v = vreset
|
||||
noise = np.sqrt(2.0*D)*np.random.randn(len(input))/np.sqrt(dt)
|
||||
for k in range(len(noise)) :
|
||||
v += (input[k]+noise[k])*dt/tau
|
||||
if v >= vthresh :
|
||||
v = vreset
|
||||
times.append(k*dt)
|
||||
spikes.append( times )
|
||||
return spikes
|
||||
|
||||
|
||||
def oupifspikes(rate, trials, duration, dt, D, drate, tau):
|
||||
# OU noise:
|
||||
rng = np.random.RandomState(54637281)
|
||||
time = np.arange(0.0, duration, dt)
|
||||
x = np.zeros(time.shape)+rate
|
||||
n = rng.randn(len(time))*drate*tau/np.sqrt(dt) + rate
|
||||
for k in range(1,len(x)) :
|
||||
x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau
|
||||
x[x<0.0] = 0.0
|
||||
spikes = pifspikes(x, trials, dt, D)
|
||||
return spikes
|
||||
|
||||
|
||||
def plot_count_fano(ax1, ax2, spikes):
|
||||
wins = np.logspace(-2, 0.0, 200)
|
||||
mean_counts = np.zeros(len(wins))
|
||||
var_counts = np.zeros(len(wins))
|
||||
for k, win in enumerate(wins):
|
||||
counts = []
|
||||
for times in spikes:
|
||||
c, _ = np.histogram(times, np.arange(0.0, duration, win))
|
||||
counts.extend(c)
|
||||
mean_counts[k] = np.mean(counts)
|
||||
var_counts[k] = np.var(counts)
|
||||
ax1.plot(mean_counts, var_counts, **lsA)
|
||||
ax1.set_xlabel('Mean count')
|
||||
ax1.set_xlim(0.0, 20.0)
|
||||
ax1.set_ylim(0.0, 20.0)
|
||||
ax2.plot(1000.0*wins, var_counts/mean_counts, **lsB)
|
||||
ax2.set_xlabel('Window', 'ms')
|
||||
ax2.set_ylim(0.0, 1.2)
|
||||
ax2.set_xscale('log')
|
||||
ax2.set_xticks([10, 100, 1000])
|
||||
ax2.set_xticklabels(['10', '100', '1000'])
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
homspikes = hompoisson(rate, trials, duration)
|
||||
inhspikes = oupifspikes(rate, trials, duration, dt, 0.3, drate, tau)
|
||||
fig, axs = plt.subplots(2, 2)
|
||||
fig.subplots_adjust(**adjust_fs(fig, top=0.5, right=1.5))
|
||||
plot_count_fano(axs[0,0], axs[0,1], homspikes)
|
||||
axs[0,0].text(0.1, 0.95, 'Poisson', transform=axs[0,0].transAxes)
|
||||
axs[0,0].set_xlabel('')
|
||||
axs[0,1].set_xlabel('')
|
||||
axs[0,0].xaxis.set_major_formatter(mpt.NullFormatter())
|
||||
axs[0,1].xaxis.set_major_formatter(mpt.NullFormatter())
|
||||
plot_count_fano(axs[1,0], axs[1,1], inhspikes)
|
||||
axs[1,0].text(0.1, 0.95, 'OU noise', transform=axs[1,0].transAxes)
|
||||
fig.text(0.01, 0.58, 'Count variance', va='center', rotation='vertical')
|
||||
fig.text(0.53, 0.58, 'Fano factor', va='center', rotation='vertical')
|
||||
plt.savefig('fanoexamples.pdf')
|
||||
plt.close()
|
||||
@ -187,40 +187,79 @@ with itself and is always 1.
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
\section{Count statistics}
|
||||
|
||||
% \begin{figure}[t]
|
||||
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill
|
||||
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
|
||||
% \titlecaption{\label{countstatsfig}Count statistic.}{}
|
||||
% \end{figure}
|
||||
Counting the number of events $n_i$ (counts) in time windows $i$ of duration $W$
|
||||
yields positive integer random numbers that are commonly quantified
|
||||
using the following measures:
|
||||
\begin{itemize}
|
||||
\item Histogram of the counts $n_i$.
|
||||
\item Average number of counts: $\mu_n = \langle n \rangle$.
|
||||
\item Variance of counts: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
|
||||
\item \entermde{Fano Faktor}{Fano factor} (variance of counts divided by average count): $F = \frac{\sigma_n^2}{\mu_n}$.
|
||||
\end{itemize}
|
||||
Of particular interest is the \enterm[firing rate!average]{average
|
||||
firing rate} $r$ (spike count per time interval, \determ{Feuerrate})
|
||||
that is given in Hertz
|
||||
\begin{figure}[t]
|
||||
\includegraphics{countexamples}
|
||||
\titlecaption{\label{countstatsfig}Count statistics.}{The
|
||||
distribution of the number of events $k$ (blue) counted within
|
||||
windows of 20\,ms (left) or 200\,ms duration (right) of the
|
||||
homogeneous Poisson spike train with a rate of 20\,Hz shown in
|
||||
\figref{rasterexamplesfig}. For Poisson spike trains these
|
||||
distributions are given by Poisson distributions (red).}
|
||||
\end{figure}
|
||||
|
||||
The most commonly used measure for characterizing spike trains is the
|
||||
\enterm[firing rate!average]{average firing rate}. The firing rate $r$
|
||||
is the average number of spikes counted within some time interval $W$
|
||||
\begin{equation}
|
||||
\label{firingrate}
|
||||
r = \frac{\langle n \rangle}{W} \; .
|
||||
r = \frac{\langle n \rangle}{W}
|
||||
\end{equation}
|
||||
and is neasured in Hertz. The average of the spike counts is taken
|
||||
over trials. For stationary spike trains (no change in statistics, in
|
||||
particular the firing rate, over time), the firing rate based on the
|
||||
spike count equals the inverse average interspike interval
|
||||
$1/\mu_{ISI}$.
|
||||
|
||||
The firing rate based on an averaged spike counts is one example of
|
||||
many statistics based on event counts. Stationary spike trains can be
|
||||
split into many segments $i$, each of duration $W$, and the number of
|
||||
events $n_i$ in each of the segments can be counted. The integer event
|
||||
counts can be quantified in the usual ways:
|
||||
\begin{itemize}
|
||||
\item Histogram of the counts $n_i$ (\figref{countstatsfig}).
|
||||
\item Average number of counts: $\mu_n = \langle n \rangle$.
|
||||
\item Variance of counts:
|
||||
$\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
|
||||
\end{itemize}
|
||||
Because spike counts are unitless and positive numbers, the
|
||||
\begin{itemize}
|
||||
\item \entermde{Fano Faktor}{Fano factor} (variance of counts divided
|
||||
by average count): $F = \frac{\sigma_n^2}{\mu_n}$.
|
||||
\end{itemize}
|
||||
is an additional measure quantifying event counts.
|
||||
|
||||
Note that all of these statistics depend on the chosen window length
|
||||
$W$. The average spike count, for example, grows linearly with $W$ for
|
||||
sufficiently large time windows: $\langle n \rangle = r W$,
|
||||
\eqnref{firingrate}. Doubling the counting window doubles the spike
|
||||
count. As does the spike-count variance (\figref{fanofig}). At smaller
|
||||
time windows the statistics of the event counts might depend on the
|
||||
particular duration of the counting window. There might be an optimal
|
||||
time window for which the variance of the spike count is minimal. The
|
||||
Fano factor plotted as a function of the time window illustrates such
|
||||
properties of point processes (\figref{fanofig}).
|
||||
|
||||
This also has consequences for information transmission in neural
|
||||
systems. The lower the variance in spike count relative to the
|
||||
averaged count, the higher the signal-to-noise ratio at which
|
||||
information encoded in the mean spike count is transmitted.
|
||||
|
||||
% \begin{figure}[t]
|
||||
% \begin{minipage}[t]{0.49\textwidth}
|
||||
% Poisson process $\lambda=100$\,Hz:\\
|
||||
% \includegraphics[width=1\textwidth]{poissonfano100hz}
|
||||
% \end{minipage}
|
||||
% \hfill
|
||||
% \begin{minipage}[t]{0.49\textwidth}
|
||||
% LIF $I=10$, $\tau_{adapt}=100$\,ms:\\
|
||||
% \includegraphics[width=1\textwidth]{lifadaptfano10-100ms}
|
||||
% \end{minipage}
|
||||
% \titlecaption{\label{fanofig}Fano factor.}{}
|
||||
% \end{figure}
|
||||
\begin{figure}[t]
|
||||
\includegraphics{fanoexamples}
|
||||
\titlecaption{\label{fanofig}
|
||||
Count variance and Fano factor.}{Variance of event counts as a
|
||||
function of mean counts obtained by varying the duration of the
|
||||
count window (left). Dividing the count variance by the respective
|
||||
mean results in the Fano factor that can be plotted as a function
|
||||
of the count window (right). For Poisson spike trains the variance
|
||||
always equals the mean counts and consequently the Fano factor
|
||||
equals one irrespective of the count window (top). A spike train
|
||||
with positive correlations between interspike intervals (caused by
|
||||
Ohrnstein-Uhlenbeck noise) has a minimum in the Fano factor, that
|
||||
is an analysis window for which the relative count variance is
|
||||
minimal somewhere close to the correlation time scale of the
|
||||
interspike intervals (bottom).}
|
||||
\end{figure}
|
||||
|
||||
\begin{exercise}{counthist.m}{}
|
||||
Implement a function \varcode{counthist()} that calculates and plots
|
||||
|
||||
Reference in New Issue
Block a user