From 84b14d72ace1172df5551062c93d55c46de19fe4 Mon Sep 17 00:00:00 2001 From: Jan Benda Date: Tue, 19 Jan 2021 00:10:39 +0100 Subject: [PATCH] [pointprocesses] moved exercise codes to exercise folder --- pointprocesses/code/counts.m | 25 ++++++ pointprocesses/code/plotcounthist.m | 2 +- pointprocesses/code/spikecounts.m | 30 ------- pointprocesses/code/spikecountshists.m | 17 ---- pointprocesses/exercises/counts.m | 25 ++++++ pointprocesses/exercises/fanoplot.m | 32 ++++++++ pointprocesses/exercises/fanoplots.m | 9 ++ pointprocesses/exercises/isihist.m | 29 +++++++ pointprocesses/exercises/isis.m | 16 ++++ pointprocesses/exercises/isiserialcorr.m | 21 +++++ pointprocesses/exercises/plotisihist.m | 31 +++++++ .../{code => exercises}/plotisihs.m | 0 pointprocesses/exercises/plotisiserialcorr.m | 14 ++++ pointprocesses/exercises/plotpoissonisih.m | 12 +++ .../{code => exercises}/plotserialcorr.m | 8 +- .../{code => exercises}/plotspikeraster.m | 0 pointprocesses/exercises/pointprocesses-1.tex | 82 +++++++++++++------ pointprocesses/exercises/pointprocesses-2.tex | 4 +- pointprocesses/exercises/rasterplot.m | 43 ++++++++++ pointprocesses/exercises/spikecountshists.m | 22 +++++ 20 files changed, 341 insertions(+), 81 deletions(-) create mode 100644 pointprocesses/code/counts.m delete mode 100644 pointprocesses/code/spikecounts.m delete mode 100644 pointprocesses/code/spikecountshists.m create mode 100644 pointprocesses/exercises/counts.m create mode 100644 pointprocesses/exercises/fanoplot.m create mode 100644 pointprocesses/exercises/fanoplots.m create mode 100644 pointprocesses/exercises/isihist.m create mode 100644 pointprocesses/exercises/isis.m create mode 100644 pointprocesses/exercises/isiserialcorr.m create mode 100644 pointprocesses/exercises/plotisihist.m rename pointprocesses/{code => exercises}/plotisihs.m (100%) create mode 100644 pointprocesses/exercises/plotisiserialcorr.m create mode 100644 pointprocesses/exercises/plotpoissonisih.m rename pointprocesses/{code => exercises}/plotserialcorr.m (53%) rename pointprocesses/{code => exercises}/plotspikeraster.m (100%) create mode 100644 pointprocesses/exercises/rasterplot.m create mode 100644 pointprocesses/exercises/spikecountshists.m diff --git a/pointprocesses/code/counts.m b/pointprocesses/code/counts.m new file mode 100644 index 0000000..90e53ef --- /dev/null +++ b/pointprocesses/code/counts.m @@ -0,0 +1,25 @@ +function n = counts(spikes, w) +% Count spikes in time windows. +% +% Arguments: +% spikes: a cell array of vectors of spike times in seconds +% w: duration of window in seconds for computing the counts +% +% Returns: +% n: vector with spike counts + tmax = spikes{1}(end); + n = []; + 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 +end diff --git a/pointprocesses/code/plotcounthist.m b/pointprocesses/code/plotcounthist.m index c487e9a..ddfa7fc 100644 --- a/pointprocesses/code/plotcounthist.m +++ b/pointprocesses/code/plotcounthist.m @@ -1,4 +1,4 @@ -function counthist(spikes, w) +function plotcounthist(spikes, w) % Plot histogram of spike counts. % % Arguments: diff --git a/pointprocesses/code/spikecounts.m b/pointprocesses/code/spikecounts.m deleted file mode 100644 index c3eef4f..0000000 --- a/pointprocesses/code/spikecounts.m +++ /dev/null @@ -1,30 +0,0 @@ -function counts = spikecounts(spikes, w) -% Compute vector of spike counts. -% -% counts = spikecounts(spikes, w) -% -% Arguments: -% spikes: a cell array of vectors of spike times in seconds -% w: observation window duration in seconds for computing the counts -% -% Returns: -% counts: vector of spike counts - - % collect spike counts: - tmax = spikes{1}(end); - counts = []; - for k = 1:length(spikes) - times = spikes{k}; -% method 1: count the number of spikes in each window: -% for tk = 0:w:tmax-w -% nn = sum((times >= tk) & (times < tk+w)); -% %nn = length(times((times >= tk) & (times < tk+w))); -% %nn = length(find((times >= tk) & (times < tk+w))); -% counts = [counts nn]; -% end -% method 2: use the hist() function to do that! - tbins = 0.5*w:w:tmax-0.5*w; - nn = hist(times, tbins); - counts = [counts nn]; - end -end diff --git a/pointprocesses/code/spikecountshists.m b/pointprocesses/code/spikecountshists.m deleted file mode 100644 index 5876386..0000000 --- a/pointprocesses/code/spikecountshists.m +++ /dev/null @@ -1,17 +0,0 @@ -w = 0.1; -bins = 0.0:1.0:10.0; - -counts{1} = spikecounts(poissonspikes, w); -counts{2} = spikecounts(pifouspikes, w); -counts{3} = spikecounts(lifadaptspikes, w); -titles = {'Poisson', 'PIF OU', 'LIF adapt'}; -for k = 1:3 - subplot(1, 3, k); - [h, b] = hist(counts{k}, bins); - bar(b, h/sum(h)); - title(titles{k}) - xlabel('Spike count'); - xlim([-0.5, 8.5]); - ylim([0.0 0.8]); -end -savefigpdf(gcf, 'spikecountshists.pdf', 20, 7); diff --git a/pointprocesses/exercises/counts.m b/pointprocesses/exercises/counts.m new file mode 100644 index 0000000..90e53ef --- /dev/null +++ b/pointprocesses/exercises/counts.m @@ -0,0 +1,25 @@ +function n = counts(spikes, w) +% Count spikes in time windows. +% +% Arguments: +% spikes: a cell array of vectors of spike times in seconds +% w: duration of window in seconds for computing the counts +% +% Returns: +% n: vector with spike counts + tmax = spikes{1}(end); + n = []; + 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 +end diff --git a/pointprocesses/exercises/fanoplot.m b/pointprocesses/exercises/fanoplot.m new file mode 100644 index 0000000..671ced9 --- /dev/null +++ b/pointprocesses/exercises/fanoplot.m @@ -0,0 +1,32 @@ +function fanoplot(spikes, titles) +% Plots fano factor as a function of window size. +% +% Arguments: +% 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 = zeros(1, length(windows)); + vc = zeros(1, length(windows)); + for j = 1:length(windows) + w = windows(j); + spikecounts = counts(spikes, w); + mc(j) = mean(spikecounts); + vc(j) = var(spikecounts); + 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, vc ./ mc, '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 + diff --git a/pointprocesses/exercises/fanoplots.m b/pointprocesses/exercises/fanoplots.m new file mode 100644 index 0000000..949eb91 --- /dev/null +++ b/pointprocesses/exercises/fanoplots.m @@ -0,0 +1,9 @@ +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 diff --git a/pointprocesses/exercises/isihist.m b/pointprocesses/exercises/isihist.m new file mode 100644 index 0000000..a24989c --- /dev/null +++ b/pointprocesses/exercises/isihist.m @@ -0,0 +1,29 @@ +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 +% +% 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); + % normalization (integral = 1): + pdf = nelements / sum(nelements) / binwidth; +end + diff --git a/pointprocesses/exercises/isis.m b/pointprocesses/exercises/isis.m new file mode 100644 index 0000000..d5eb7d6 --- /dev/null +++ b/pointprocesses/exercises/isis.m @@ -0,0 +1,16 @@ +function isivec = isis(spikes) +% returns a single list of isis computed from all trials in spikes +% +% Arguments: +% spikes: a cell array of vectors of spike times in seconds +% +% Returns: +% isivec: a column vector with all the interspike intervals + isivec = []; + for k = 1:length(spikes) + difftimes = diff(spikes{k}); + % difftimes(:) ensures a column vector + % regardless of the type of vector in spikes{k} + isivec = [isivec; difftimes(:)]; + end +end diff --git a/pointprocesses/exercises/isiserialcorr.m b/pointprocesses/exercises/isiserialcorr.m new file mode 100644 index 0000000..fc5ec17 --- /dev/null +++ b/pointprocesses/exercises/isiserialcorr.m @@ -0,0 +1,21 @@ +function [isicorr, lags] = isiserialcorr(isivec, maxlag) +% serial correlation of interspike intervals +% +% Arguments: +% isivec: vector of interspike intervals in seconds +% maxlag: the maximum lag +% +% Returns: +% isicorr: vector with the serial correlations for lag 0 to maxlag +% lags: vector with the lags corresponding to isicorr + lags = 0:maxlag; + isicorr = zeros(size(lags)); + for k = 1:length(lags) + lag = lags(k); + if length(isivec) > lag+10 % ensure "enough" data + % NOTE: the arguments to corr must be column vectors! + % We insure this already in the isis() function. + isicorr(k) = corr(isivec(1:end-lag), isivec(lag+1:end)); + end + end +end diff --git a/pointprocesses/exercises/plotisihist.m b/pointprocesses/exercises/plotisihist.m new file mode 100644 index 0000000..2f3cf93 --- /dev/null +++ b/pointprocesses/exercises/plotisihist.m @@ -0,0 +1,31 @@ +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: + 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), ... + 'Units', 'normalized', 'HorizontalAlignment', 'right') + text(0.95, 0.7, sprintf('std=%.1f ms', 1000.0*sdisi), ... + 'Units', 'normalized', 'HorizontalAlignment', 'right') + text(0.95, 0.6, sprintf('CV=%.2f', sdisi/misi), ... + 'Units', 'normalized', 'HorizontalAlignment', 'right') +end + diff --git a/pointprocesses/code/plotisihs.m b/pointprocesses/exercises/plotisihs.m similarity index 100% rename from pointprocesses/code/plotisihs.m rename to pointprocesses/exercises/plotisihs.m diff --git a/pointprocesses/exercises/plotisiserialcorr.m b/pointprocesses/exercises/plotisiserialcorr.m new file mode 100644 index 0000000..ce806e3 --- /dev/null +++ b/pointprocesses/exercises/plotisiserialcorr.m @@ -0,0 +1,14 @@ +function isicorr = plotisiserialcorr(isivec, maxlag) +% plot serial correlation of interspike intervals +% +% Arguments: +% isivec: vector of interspike intervals in seconds +% maxlag: the maximum lag + [isicorr, lags] = isiserialcorr(isivec, maxlag); + plot(lags, isicorr, '-b'); + hold on; + scatter(lags, isicorr, 20.0, 'b', 'filled'); + hold off; + xlabel('Lag k') + ylabel('\rho_k') +end diff --git a/pointprocesses/exercises/plotpoissonisih.m b/pointprocesses/exercises/plotpoissonisih.m new file mode 100644 index 0000000..2218196 --- /dev/null +++ b/pointprocesses/exercises/plotpoissonisih.m @@ -0,0 +1,12 @@ +maxisi = 300.0; +poissonisis = isis(poissonspikes); +rate = 1.0/mean(poissonisis); +plotisihist(poissonisis, 0.001); +tt = linspace(0.0, 0.001*maxisi, 200); +pexp = rate*exp(-tt*rate); +hold on; +plot(1000.0*tt, pexp, 'r') +hold off; +xlim([0, maxisi]) +title('Poisson'); +savefigpdf(gcf, 'poissonisihist.pdf', 10, 7); diff --git a/pointprocesses/code/plotserialcorr.m b/pointprocesses/exercises/plotserialcorr.m similarity index 53% rename from pointprocesses/code/plotserialcorr.m rename to pointprocesses/exercises/plotserialcorr.m index 05f3ea8..ab33e50 100644 --- a/pointprocesses/code/plotserialcorr.m +++ b/pointprocesses/exercises/plotserialcorr.m @@ -1,17 +1,17 @@ maxlag = 10; rrange = [-0.5, 1.05]; subplot(1, 3, 1); -isiserialcorr(poissonisis, maxlag); +plotisiserialcorr(poissonisis, maxlag); ylim(rrange) title('Poisson'); subplot(1, 3, 2); -isiserialcorr(pifouisis, maxlag); +plotisiserialcorr(pifouisis, maxlag); ylim(rrange) title('PIF OU'); subplot(1, 3, 3); -isiserialcorr(lifadaptisis, maxlag); +plotisiserialcorr(lifadaptisis, maxlag); ylim(rrange) title('LIF adapt'); -savefigpdf(gcf, 'serialcorr.pdf', 20, 7); \ No newline at end of file +savefigpdf(gcf, 'serialcorr.pdf', 20, 7); diff --git a/pointprocesses/code/plotspikeraster.m b/pointprocesses/exercises/plotspikeraster.m similarity index 100% rename from pointprocesses/code/plotspikeraster.m rename to pointprocesses/exercises/plotspikeraster.m diff --git a/pointprocesses/exercises/pointprocesses-1.tex b/pointprocesses/exercises/pointprocesses-1.tex index 2da20f9..6edb623 100644 --- a/pointprocesses/exercises/pointprocesses-1.tex +++ b/pointprocesses/exercises/pointprocesses-1.tex @@ -15,6 +15,9 @@ \begin{questions} + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + \question \qt{Read in chapter 10 on ``Spike-train analysis'' sections 10.1 -- 10.5!}\vspace{-3ex} + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \question \qt{Statistics of interspike intervals} Download the files \code{poisson.mat}, @@ -61,11 +64,12 @@ of the spike. When appropriate, the function should use milliseconds for the time axis instead of seconds. - Use this function to plot the first second of the - spike rasters of the three neurons. + Use this function to plot the first second of the spike rasters of + the three neurons. Carefully look at the raster plots. How do the + spike rasters differ between the three neurons? \begin{solution} - \lstinputlisting{../code/rasterplot.m} - \lstinputlisting{../code/plotspikeraster.m} + \lstinputlisting{rasterplot.m} + \lstinputlisting{plotspikeraster.m} \mbox{}\\[-3ex] \colorbox{white}{\includegraphics[width=1\textwidth]{spikeraster}} \end{solution} @@ -73,7 +77,7 @@ \part Write a function that returns a single vector containing the interspike intervals of all trials of spike times. \begin{solution} - \lstinputlisting{../code/isis.m} + \lstinputlisting{isis.m} \end{solution} \part Write a function that computes an estimate of the @@ -90,33 +94,51 @@ the standard deviation and the coefficient of variation and display them in the plot as well. - Use this and the previous functions to compare the - interspike interval statistics of the three neurons. + Use this and the previous functions to compare the interspike + interval statistics of the three neurons. How do the ISI + histograms differ between the three neurons? \begin{solution} - \lstinputlisting{../code/isihist.m} - \lstinputlisting{../code/plotisihist.m} - \lstinputlisting{../code/plotisihs.m} + \lstinputlisting{isihist.m} + \lstinputlisting{plotisihist.m} + \lstinputlisting{plotisihs.m} \mbox{}\\[-3ex] \colorbox{white}{\includegraphics[width=1\textwidth]{isihist}} \end{solution} + + \part Compare the interspike interval histogram of the Poisson + spike train to the expected exponential distribution + \[ p(T) = \lambda e^{-\lambda T} \; .\] of a Poisson process with + firing rate $\lambda$. Estimate the firing rate from the inverse + of the mean interspike interval. + \begin{solution} + \lstinputlisting{plotpoissonisih.m} + %\mbox{}\\[-3ex] + %\colorbox{white}{\includegraphics[width=1\textwidth]{poissonisihist}} + \end{solution} % XXX Add return map!!! XXX - \part Write a function that computes and plots the serial - correlations of interspike intervals for lags upto - \code{maxlag}. The serial correlations $\rho_k$ for lag $k$ of the - interspike intervals $T_i$ are the correlation coefficients - between interspike intervals $T_i$ and the intervals $T_{i+k}$ - that are shifted by lag $k$: + \part Write a function that computes the serial correlations of + interspike intervals for lags up to \code{maxlag}. The serial + correlations $\rho_k$ for lag $k$ of the interspike intervals + $T_i$ are the correlation coefficients between interspike + intervals $T_i$ and the intervals $T_{i+k}$ that are shifted by + lag $k$: \[ \rho_k = \frac{\langle (T_{i+k} - \langle T \rangle)(T_i - \langle T \rangle) \rangle}{\langle (T_i - \langle T \rangle)^2\rangle} = \frac{{\rm cov}(T_{i+k}, T_i)}{{\rm - var}(T_i)} = {\rm corr}(T_{i+k}, T_i) \] + var}(T_i)} = {\rm corr}(T_{i+k}, T_i) \] + + Write another function that plots the serial correlations. - Use this function to compare the serial correlations of the - interspike intervals of the three neurons. + Use the functions to compare the serial correlations of the + interspike intervals of the three neurons. What do the serial + correlations mean? How do they relate to what you see in the + raster plots? Do the serial correlations of the Poisson spike + train match the expectations? \begin{solution} - \lstinputlisting{../code/isiserialcorr.m} - \lstinputlisting{../code/plotserialcorr.m} + \lstinputlisting{isiserialcorr.m} + \lstinputlisting{plotisiserialcorr.m} + \lstinputlisting{plotserialcorr.m} \colorbox{white}{\includegraphics[width=1\textwidth]{serialcorr}} \end{solution} @@ -127,16 +149,22 @@ \question \qt{Statistics of spike counts} Now let's have a look at the statistics of the spike counts. \begin{parts} - \part Write a function that counts and returns the number of - spikes in windows of a given width $W$. + \part Write a function that counts and returns a vector with the + number of spikes in windows of a given width $W$. Use this function to generate a properly normalized histogram of spike counts for the data of the three types of neurons. Use 100\,ms for the window width. + + Compare the distributions with the Poisson distribution expected + for a Poisson spike train: + \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \; , \] where + $\lambda$ is the rate of the spike train that you should estimate + from the data. \begin{solution} - \lstinputlisting{../code/spikecounts.m} + \lstinputlisting{counts.m} \newsolutionpage - \lstinputlisting{../code/spikecountshists.m} + \lstinputlisting{spikecounthists.m} \colorbox{white}{\includegraphics[width=1\textwidth]{spikecountshists}} \end{solution} @@ -148,9 +176,9 @@ (the mean spike count increases for larger window widths). The other plot showing the Fano factor as a function of window width. \begin{solution} - \lstinputlisting{../code/fanoplot.m} + \lstinputlisting{fanoplot.m} \newsolutionpage - \lstinputlisting{../code/fanoplots.m} + \lstinputlisting{fanoplots.m} \colorbox{white}{\includegraphics[width=1\textwidth]{fanoplotspoisson}} \colorbox{white}{\includegraphics[width=1\textwidth]{fanoplotspifou}} \colorbox{white}{\includegraphics[width=1\textwidth]{fanoplotslifadapt}} diff --git a/pointprocesses/exercises/pointprocesses-2.tex b/pointprocesses/exercises/pointprocesses-2.tex index 883b0b2..cd7a2f3 100644 --- a/pointprocesses/exercises/pointprocesses-2.tex +++ b/pointprocesses/exercises/pointprocesses-2.tex @@ -103,7 +103,7 @@ f\"ur gen\"ugend kleine $\Delta t$. soll. Hinweis: es gibt eine \code{matlab} Funktion, die die Fakult\"at $n!$ berechnet. \begin{solution} - \lstinputlisting{../code/counthist.m} + \lstinputlisting{counthist.m} \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}} \colorbox{white}{\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}} \end{solution} @@ -115,7 +115,7 @@ f\"ur gen\"ugend kleine $\Delta t$. den Ergebnissen angefertigt werden: (i) Varianz gegen Mittelwert der counts. (ii) Fano Faktor als Funktion der Fensterbreite. \begin{solution} - \lstinputlisting{../code/fano.m} + \lstinputlisting{fanoplot.m} \colorbox{white}{\includegraphics[width=0.98\textwidth]{poissonfano100hz}} \end{solution} diff --git a/pointprocesses/exercises/rasterplot.m b/pointprocesses/exercises/rasterplot.m new file mode 100644 index 0000000..ccbfe60 --- /dev/null +++ b/pointprocesses/exercises/rasterplot.m @@ -0,0 +1,43 @@ +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 + in_msecs = tmax < 1.5 + spiketimes = []; + trials = []; + ntrials = length(spikes); + for k = 1:ntrials + times = spikes{k}; + times = times(times