diff --git a/pointprocesses/lecture/countexamples.py b/pointprocesses/lecture/countexamples.py new file mode 100644 index 0000000..7e0f98f --- /dev/null +++ b/pointprocesses/lecture/countexamples.py @@ -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() diff --git a/pointprocesses/lecture/fanoexamples.py b/pointprocesses/lecture/fanoexamples.py new file mode 100644 index 0000000..dc68c0e --- /dev/null +++ b/pointprocesses/lecture/fanoexamples.py @@ -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= 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() diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex index f2286e2..d13d0ab 100644 --- a/pointprocesses/lecture/pointprocesses.tex +++ b/pointprocesses/lecture/pointprocesses.tex @@ -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