teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

[pointprocesses] improved spike count and fano factor exercises

Browse Source
This commit is contained in:
Jan Benda 2021-01-18 23:03:30 +01:00
parent 753d756409
commit 7ec4cddd62
14 changed files with 109 additions and 220 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
pointprocesses/code/counthist.m
View File

@ -1,37 +0,0 @@
function counthist(spikes, w)
% Plot histogram of spike counts.
%
% counthist(spikes, w)
%
% Arguments:
% spikes: a cell array of vectors of spike times in seconds
% w: duration of window in seconds for computing the counts
% collect spike counts:
tmax = spikes{1}(end);
n = [];
r = [];
for k = 1:length(spikes)
times = spikes{k};
% alternative 1: count the number of spikes in each window:
% for tk = 0:w:tmax-w
% nn = sum((times >= tk) & (times < tk+w));
% %nn = length(find((times >= tk) & (times < tk+w)));
% n = [n, nn];
% end
% alternative 2: use the hist() function to do that!
tbins = 0.5*w:w:tmax-0.5*w;
nn = hist(times, tbins);
n = [n, nn];
end
% histogram of spike counts:
maxn = max(n);
[counts, bins] = hist(n, 0:1:maxn+10);
% normalize to probabilities:
counts = counts / sum(counts);
% plot:
bar(bins, counts);
xlabel('counts k');
ylabel('P(k)');
end

39
pointprocesses/code/fano.m
View File

@ -1,39 +0,0 @@
function fano( spikes )
% computes fano factor as a function of window size
% spikes: a cell array of vectors of spike times
tmax = spikes{1}(end);
windows = 0.01:0.05:0.01*tmax;
mc = windows;
vc = windows;
ff = windows;
fs = windows;
for j = 1:length(windows)
w = windows( j );
% collect counts:
n = [];
for k = 1:length(spikes)
for tk = 0:w:tmax-w
nn = sum( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) );
%nn = length( find( ( spikes{k} >= tk ) & ( spikes{k} < tk+w ) ) );
n = [ n nn ];
end
end
% statistics for current window:
mc(j) = mean( n );
vc(j) = var( n );
ff(j) = vc( j )/mc( j );
fs(j) = sqrt(vc( j )/mc( j ));
end
subplot( 1, 2, 1 );
scatter( mc, vc, 'filled' );
xlabel( 'Mean count' );
ylabel( 'Count variance' );
subplot( 1, 2, 2 );
scatter( 1000.0*windows, fs, 'filled' );
xlabel( 'Window W [ms]' );
ylabel( 'Fano factor' );
end

33
pointprocesses/code/fanoplot.m
View File

@ -1,33 +0,0 @@
function fanoplot(spikes, titles)
% computes and plots fano factor as a function of window size
% spikes: a cell array of vectors of spike times
% titles: string that is used as a title for the plots
windows = logspace(-3.0, -0.5, 100);
mc = windows;
vc = windows;
ff = windows;
for j = 1:length(windows)
w = windows(j);
counts = spikecounts(spikes, w);
% statistics for current window:
mc(j) = mean(counts);
vc(j) = var(counts);
ff(j) = vc(j)/mc(j);
end
subplot(1, 2, 1);
scatter(mc, vc, 'filled');
title(titles);
xlabel('Mean count');
ylabel('Count variance');
subplot(1, 2, 2);
scatter(1000.0*windows, ff, 'filled');
title(titles);
xlabel('Window [ms]');
ylabel('Fano factor');
xlim(1000.0*[windows(1) windows(end)])
ylim([0.0 1.1]);
set(gca, 'XScale', 'log');
end

9
pointprocesses/code/fanoplots.m
View File

@ -1,9 +0,0 @@
spikes{1} = poissonspikes;
spikes{2} = pifouspikes;
spikes{3} = lifadaptspikes;
idents = {'poisson', 'pifou', 'lifadapt'};
for k = 1:3
figure(k)
fanoplot(spikes{k}, titles{k});
savefigpdf(gcf, sprintf('fanoplots%s.pdf', idents{k}), 20, 7);
end

3
pointprocesses/code/hompoissonspikes.m
View File

@ -2,8 +2,6 @@ function spikes = hompoissonspikes(rate, trials, tmax)
% Generate spike times of a homogeneous poisson process
% using the exponential interspike interval distribution.
%
% spikes = hompoissonspikes(rate, trials, tmax)
%
% Arguments:
% rate: the rate of the Poisson process in Hertz
% trials: number of trials that should be generated
@ -11,7 +9,6 @@ function spikes = hompoissonspikes(rate, trials, tmax)
%
% Returns:
% spikes: a cell array of vectors of spike times in seconds
spikes = cell(trials, 1);
mu = 1.0/rate;
nintervals = 2*round(tmax/mu);

14
pointprocesses/code/isihist.m
View File

@ -1,25 +1,13 @@
function [pdf, centers] = isihist(isis, binwidth)
% Compute normalized histogram of interspike intervals.
%
% [pdf, centers] = isihist(isis, binwidth)
%
% Arguments:
% isis: vector of interspike intervals in seconds
% binwidth: optional width in seconds to be used for the isi bins
% binwidth: width in seconds to be used for the ISI bins
%
% Returns:
% pdf: vector with pdf of interspike intervals in Hertz
% centers: vector with centers of interspikeintervalls in seconds
if nargin < 2
% compute good binwidth:
nperbin = 200; % average number of data points per bin
bins = length(isis)/nperbin; % number of bins
binwidth = max(isis)/bins;
if binwidth < 5e-4 % half a millisecond
binwidth = 5e-4;
end
end
bins = 0.5*binwidth:binwidth:max(isis);
% histogram data:
[nelements, centers] = hist(isis, bins);

4
pointprocesses/code/isis.m
View File

@ -1,13 +1,11 @@
function isivec = isis(spikes)
% returns a single list of isis computed from all trials in spikes
%
% isivec = isis(spikes)
%
% Arguments:
% spikes: a cell array of vectors of spike times in seconds
%
% Returns:
% isivec: a column vector with all the interspike intervalls
% isivec: a column vector with all the interspike intervals
isivec = [];
for k = 1:length(spikes)
difftimes = diff(spikes{k});

5
pointprocesses/code/isiserialcorr.m
View File

@ -1,8 +1,6 @@
function [isicorr, lags] = isiserialcorr(isivec, maxlag)
% serial correlation of interspike intervals
%
% isicorr = isiserialcorr(isivec, maxlag)
%
% Arguments:
% isivec: vector of interspike intervals in seconds
% maxlag: the maximum lag
@ -16,8 +14,7 @@ function [isicorr, lags] = isiserialcorr(isivec, maxlag)
lag = lags(k);
if length(isivec) > lag+10 % ensure "enough" data
% NOTE: the arguments to corr must be column vectors!
% We insure this in the isis() function that
% generates the isivec.
% We insure this already in the isis() function.
isicorr(k) = corr(isivec(1:end-lag), isivec(lag+1:end));
end
end

38
pointprocesses/code/plotcounthist.m
View File

@ -1,24 +1,14 @@
w = 0.1;
cmax = 8;
pmax = 0.5;
subplot(1, 3, 1);
counthist(poissonspikes, w);
xlim([0 cmax])
set(gca, 'XTick', 0:2:cmax)
ylim([0 pmax])
title('Poisson');
subplot(1, 3, 2);
counthist(pifouspikes, w);
xlim([0 cmax])
set(gca, 'XTick', 0:2:cmax)
ylim([0 pmax])
title('PIF OU');
subplot(1, 3, 3);
counthist(lifadaptspikes, w);
xlim([0 cmax])
set(gca, 'XTick', 0:2:cmax)
ylim([0 pmax])
title('LIF adapt');
savefigpdf(gcf, 'counthist.pdf', 20, 7);
function counthist(spikes, w)
% Plot histogram of spike counts.
%
% Arguments:
% spikes: a cell array of vectors of spike times in seconds
% w: duration of window in seconds for computing the counts
n = counts(spikes, w);
maxn = max(n);
[counts, bins] = hist(n, 0:1:maxn+10);
counts = counts / sum(counts);
bar(bins, counts);
xlabel('counts k');
ylabel('P(k)');
end

31
pointprocesses/code/plotfanofactor.m Normal file
View File

@ -0,0 +1,31 @@
function plotfanofactor(spikes, wmin, wmax)
% Compute and plot Fano factor as a function of window size.
%
% Arguments:
% spikes: a cell array of vectors of spike times in seconds
% wmin: minimum window size in seconds
% wmax: maximum window size in seconds
windows = logspace(log10(wmin), log10(wmax), 100);
mc = zeros(1, length(windows));
vc = zeros(1, length(windows));
for k = 1:length(windows)
w = windows(k);
n = counts(spikes, w);
mc(k) = mean(n);
vc(k) = var(n);
end
subplot(1, 2, 1);
scatter(mc, vc, 'filled');
xlabel('Mean count');
ylabel('Count variance');
subplot(1, 2, 2);
scatter(1000.0*windows, vc ./ mc, 'filled');
xlabel('Window [ms]');
ylabel('Fano factor');
xlim(1000.0*[windows(1) windows(end)])
ylim([0.0 1.1]);
set(gca, 'XScale', 'log');
end

15
pointprocesses/code/plotisihist.m
View File

@ -1,24 +1,13 @@
function plotisihist(isis, binwidth)
% Plot and annotate histogram of interspike intervals.
%
% plotisihist(isis, binwidth)
%
% Arguments:
% isis: vector of interspike intervals in seconds
% binwidth: optional width in seconds to be used for the isi bins
% compute normalized histogram:
if nargin < 2
[pdf, centers] = isihist(isis);
else
[pdf, centers] = isihist(isis, binwidth);
end
% plot:
% binwidth: width in seconds to be used for the ISI bins
[pdf, centers] = isihist(isis, binwidth);
bar(1000.0*centers, pdf); % milliseconds on x-axis
xlabel('ISI [ms]')
ylabel('p(ISI) [1/s]')
% annotation:
misi = mean(isis);
sdisi = std(isis);
text(0.95, 0.8, sprintf('mean=%.1f ms', 1000.0*misi), ...

2
pointprocesses/code/plotisiserialcorr.m
View File

@ -1,8 +1,6 @@
function isicorr = plotisiserialcorr(isivec, maxlag)
% plot serial correlation of interspike intervals
%
% plotisiserialcorr(isivec, maxlag)
%
% Arguments:
% isivec: vector of interspike intervals in seconds
% maxlag: the maximum lag

4
pointprocesses/code/rasterplot.m
View File

@ -1,8 +1,6 @@
function rasterplot(spikes, tmax)
% Display a spike raster of the spike times given in spikes.
%
% rasterplot(spikes, tmax)
%
% Arguments:
% spikes: a cell array of vectors of spike times in seconds
% tmax: plot spike raster up to tmax seconds
@ -21,7 +19,7 @@ function rasterplot(spikes, tmax)
ones(1, length(times)) * (k+0.4); ...
ones(1, length(times)) * nan]];
end
% convert matrices into simple vectors of (x,y) pairs:
% convert matrices into column vectors of (x,y) pairs:
spiketimes = spiketimes(:);
trials = trials(:);
% plotting this is lightning fast:

95
pointprocesses/lecture/pointprocesses.tex
View File

@ -239,19 +239,18 @@ describing univariate data sets of real numbers:
\begin{exercise}{isihist.m}{}
Implement a function \varcode{isihist()} that calculates the
normalized interspike interval histogram. The function should take
two input arguments; (i) a vector of interspike intervals and (ii)
the width of the bins used for the histogram. It returns the
normalized interspike interval histogram. The function should take a
vector of interspike intervals and the width of the bins to be used
for the histogram as input arguments. The function returns the
probability density as well as the centers of the bins.
\end{exercise}
\begin{exercise}{plotisihist.m}{}
Implement a function that takes the returned values of
\varcode{isihist()} as input arguments and then plots the data. The
plot should show the histogram with the x-axis scaled to
milliseconds and should be annotated with the average ISI, the
standard deviation, and the coefficient of variation of the ISIs
(\figref{isihexamplesfig}).
Implement a function that uses the \varcode{isihist()} function from
the previous exercise to plot an ISI histogram. The plot shows the
histogram with the x-axis scaled to milliseconds, annotated with the
average ISI, the standard deviation, and the coefficient of
variation of the ISIs (\figref{isihexamplesfig}).
\end{exercise}
\subsection{Interval correlations}
@ -332,8 +331,8 @@ characterized by the non-zero serial correlations
\begin{exercise}{isiserialcorr.m}{}
Implement a function \varcode{isiserialcorr()} that takes a vector
of interspike intervals as input argument and calculates the serial
correlations up to some maximum lag.
of interspike intervals as input and computes serial correlations up
to some maximum lag.
\end{exercise}
\begin{exercise}{plotisiserialcorr.m}{}
@ -398,27 +397,10 @@ Because spike counts are unitless and positive numbers the
\end{equation}
is a commonly used measure for quantifying the variability of event
counts relative to the mean number of events. In particular for
homogeneous Poisson processes the Fano factor equals one,
independently of the time window $W$.
homogeneous Poisson processes the Fano factor equals exactly one and
is independent of the time window $W$.
\end{itemize}
Note that all of these statistics depend in general 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]
\includegraphics{fanoexamples}
\titlecaption{\label{fanofig} Fano factor.}{Counting events in time
@ -435,13 +417,52 @@ information encoded in the mean spike count is transmitted.
interspike intervals (right).}
\end{figure}
\begin{exercise}{counthist.m}{}
Implement a function \varcode{counthist()} that calculates and plots
the distribution of spike counts observed in a certain time
window. The function should take two input arguments: a cell-array
of vectors containing the spike times in seconds observed in a
number of trials, and the duration of the time window that is used
to evaluate the counts.
Note that all of these statistics depend in general on the chosen
window length $W$ used for counting the events. 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 could 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 in a single graph (\figref{fanofig}).
This also has consequences for information transmission in neural
systems. The membrane time constant of a post-synaptic neuron defines
a counting window. If input spikes arrive within about one time
constant they are added up. 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. If the membrane time constant of the target neuron
matches the counting window where the Fano factor is minimal, then
information is potentially transmitted at the highest possible
signal-to-noise ratio. For this reason, the Fano factor is used in the
Neurosciences to quantify and analyze reliability of neuronal
responses.
\begin{exercise}{counts.m}{}
Write a function \varcode{counts()} that counts the number of spikes
in windows of given duration and returns the counts in a single
vector. Spike times are passed as a cell-array of vectors,
containing the spike times in seconds observed in a number of
trials, to the function.
\end{exercise}
\begin{exercise}{plotcounthist.m}{}
Write a function that takes a cell-array with spike times as input
and plots a normalized histogram of the spike counts counted in
windows of a given duration.
\end{exercise}
\begin{exercise}{plotfanofactor.m}{}
Write a function that takes a cell-array with spike times as input
and plots in one plot count variances a function of the
corresponding mean counts and in a second plot the Fano factor as a
function of the duration of the count window in logarithmic
scale. Two arguments of the function take the minimum and maximum
duration of the count window.
\end{exercise}