Merge branch 'master' of raven.am28.uni-tuebingen.de:scientificComputing

Conflicts:
	pointprocesses/code/binnedRate.m
	pointprocesses/code/convolutionRate.m
	pointprocesses/code/instantaneousRate.m
	programming/lecture/programming.tex
	programmingstyle/code/calculateSines.m
	programmingstyle/lecture/programmingstyle.tex
This commit is contained in:
2015-11-30 12:03:30 +01:00
66 changed files with 1736 additions and 1549 deletions

View File

@@ -1,19 +1,24 @@
function [time, rate] = binnedRate(spike_times, bin_width, dt, t_max)
% Calculates the firing rate with the binning method. The hist
% function is used to count the number of spikes in each bin.
% Arguments:
% spike_times, vector containing the times of the spikes.
% bin_width, the width of the bins in seconds.
% dt, the temporal resolution.
% t_max, the tiral duration.
function [time, rate] = binned_rate(spikes, bin_width, dt, t_max)
% PSTH computed with binning method.
% The hist funciton is used to count the number of spikes in each bin.
%
% Returns two vectors containing the time and the rate.
% [time, rate] = binned_rate(spikes, bin_width, dt, t_max)
%
% Arguments:
% spikes : vector containing the times of the spikes.
% bin_width: the width of the bins in seconds.
% dt : the temporal resolution.
% t_max : the tiral duration.
%
% Returns:
% two vectors containing the time and the rate.
time = 0:dt:t_max-dt;
bins = 0:bin_width:t_max;
rate = zeros(size(time));
h = hist(spike_times, bins) ./ bin_width;
for i = 2:length(bins)
rate(round(bins(i - 1) / dt) + 1:round(bins(i) / dt)) = h(i);
time = 0:dt:t_max-dt;
bins = 0:bin_width:t_max;
rate = zeros(size(time));
h = hist(spikes, bins) ./ bin_width;
for i = 2:length(bins)
rate(round(bins(i - 1) / dt) + 1:round(bins(i) / dt)) = h(i);
end
end

View File

@@ -1,17 +1,20 @@
function [time, rate] = convolutionRate(spike_times, sigma, dt, t_max)
% Calculates the firing rate with the convolution method.
% Arguments:
% spike_times, a vector containing the spike times.
% sigma, the standard deviation of the Gaussian kernel
% in seconds.
% dt, the temporal resolution in seconds.
% t_max, the trial duration in seconds.
function [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
% PSTH computed with convolution method.
%
% Returns two vectors containing the time and the rate.
% [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
%
% Arguments:
% spikes: a vector containing the spike times.
% sigma : the standard deviation of the Gaussian kernel in seconds.
% dt : the temporal resolution in seconds.
% t_max : the trial duration in seconds.
%
% Returns:
two vectors containing the time and the rate.
time = 0:dt:t_max - dt;
rate = zeros(size(time));
spike_indices = round(spike_times / dt);
spike_indices = round(spikes / dt);
rate(spike_indices) = 1;
kernel = gaussKernel(sigma, dt);

View File

@@ -1,22 +1,24 @@
function [time, rate] = instantaneousRate(spike_times, dt, t_max)
% Function calculates the firing rate as the inverse of the interspike
% interval.
% Arguments:
% spike_times, vector containing the times of the spikes.
% dt, the temporal resolutions of the recording.
% t_max, the duration of the trial.
function [time, rate] = instantaneous_rate(spikes, dt, t_max)
% Firing rate as the inverse of the interspike interval.
%
% Returns two vectors containing the time and the rate.
% [time, rate] = instantaneous_rate(spikes, dt, t_max)
%
% Arguments:
% spikes: vector containing the times of the spikes.
% dt : the temporal resolutions of the recording.
% t_max : the duration of the trial.
%
% Returns:
% the vector representing time and a vector containing the rate.
time = 0:dt:t_max-dt;
rate = zeros(size(time));
time = 0:dt:t_max-dt;
rate = zeros(size(time));
isis = diff([0 spike_times]);
inst_rate = 1 ./ isis;
spike_indices = [1 round(spike_times ./ dt)];
isis = diff([0 spikes]);
inst_rate = 1 ./ isis;
spike_indices = [1 round(spikes ./ dt)];
for i = 2:length(spike_indices)
rate(spike_indices(i - 1):spike_indices(i)) = inst_rate(i - 1);
for i = 2:length(spike_indices)
rate(spike_indices(i - 1):spike_indices(i)) = inst_rate(i - 1);
end
end

View File

@@ -1,7 +1,7 @@
function [pdf, centers] = isi_hist(isis, binwidth)
function [pdf, centers] = isiHist(isis, binwidth)
% Compute normalized histogram of interspike intervals.
%
% [pdf, centers] = isi_hist(isis, binwidth)
% [pdf, centers] = isiHist(isis, binwidth)
%
% Arguments:
% isis: vector of interspike intervals in seconds

View File

@@ -2,7 +2,12 @@ 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
% isivec: a column vector with all the interspike intervalls
%
% Returns:
% isivec: a column vector with all the interspike intervalls
isivec = [];
@@ -13,4 +18,3 @@ function isivec = isis( spikes )
isivec = [ isivec; difftimes(:) ];
end
end

View File

@@ -1,7 +1,7 @@
function plot_isi_hist(isis, binwidth)
function plotISIHist(isis, binwidth)
% Plot and annotate histogram of interspike intervals.
%
% isihist(isis, binwidth)
% plotISIHist(isis, binwidth)
%
% Arguments:
% isis: vector of interspike intervals in seconds
@@ -9,9 +9,9 @@ function plot_isi_hist(isis, binwidth)
% compute normalized histogram:
if nargin < 2
[pdf, centers] = isi_hist(isis);
[pdf, centers] = isiHist(isis);
else
[pdf, centers] = isi_hist(isis, binwidth);
[pdf, centers] = isiHist(isis, binwidth);
end
% plot:

View File

@@ -1,19 +1,19 @@
function s_est = reconstructStimulus(spike_times, sta, stim_duration, dt)
% Function estimates the stimulus from the Spike-Triggered-Average
% (sta).
function s_est = reconstructStimulus(spikes, sta, duration, deltat)
% Estimate the stimulus from the spike-triggered-average (STA).
%
% s_est = reconstructStimulus(spikes, sta, duration, deltat)
%
% Arguments:
% spike_times, a vector containing the spike times in seconds.
% sta, a vector containing the spike-triggered-average.
% stim_duration, the total duration of the stimulus.
% dt, the sampling interval given in seconds.
% spikes : a vector containing the spike times in seconds.
% sta : a vector containing the spike-triggered-average.
% duration: the total duration of the stimulus.
% deltat : the time step of the stimulus in seconds.
%
% Returns:
% the estimated stimulus.
s_est = zeros(round(stim_duration / dt), 1);
binary_spikes = zeros(size(s_est));
binary_spikes(round(spike_times ./ dt)) = 1;
s_est = conv(binary_spikes, sta, 'same');
% s_est: vector with the estimated stimulus.
s_est = zeros(round(duration / deltat), 1);
binary_spikes = zeros(size(s_est));
binary_spikes(round(spikes ./ deltat)) = 1;
s_est = conv(binary_spikes, sta, 'same');
end

View File

@@ -1,32 +1,31 @@
function [sta, std_sta, valid_spikes] = spikeTriggeredAverage(stimulus, spike_times, count, sampling_rate)
% Function estimates the Spike-Triggered-Average (sta).
function [sta, std_sta, n_spikes] = spikeTriggeredAverage(stimulus, spikes, count, deltat)
% Estimate the spike-triggered-average (STA).
%
% [sta, std_sta, n_spikes] = spikeTriggeredAverage(stimulus, spikes, count, deltat)
%
% Arguments:
% stimulus, a vector containing stimulus intensities
% as a function of time.
% spike_times, a vector containing the spike times
% in seconds.
% count, the number of datapoints that are taken around
% the spike times.
% sampling_rate, the sampling rate of the stimulus.
% stimulus: vector of stimulus intensities as a function of time.
% spikes : vector with spike times in seconds.
% count : number of datapoints that are taken around the spike times.
% deltat : the time step of the stimulus in seconds.
%
% Returns:
% the sta, a vector containing the staandard deviation and
% the number of spikes taken into account.
% sta : vector with the STA.
% std_sta : standard deviation of the STA.
% n_spikes: number of spikes contained in STA.
snippets = zeros(numel(spike_times), 2*count);
valid_spikes = 1;
for i = 1:numel(spike_times)
t = spike_times(i);
index = round(t*sampling_rate);
if index <= count || (index + count) > length(stimulus)
continue
snippets = zeros(numel(spikes), 2*count);
n_spikes = 0;
for i = 1:numel(spikes)
t = spikes(i);
index = round(t/deltat);
if index <= count || (index + count) > length(stimulus)
continue
end
snippets(n_spikes,:) = stimulus(index-count:index+count-1);
n_spikes = n_spikes + 1;
end
snippets(valid_spikes,:) = stimulus(index-count:index+count-1);
valid_spikes = valid_spikes + 1;
snippets(n_spikes+1:end,:) = [];
sta = mean(snippets, 1);
std_sta = std(snippets,[],1);
end
snippets(valid_spikes:end,:) = [];
sta = mean(snippets, 1);
std_sta = std(snippets,[],1);

View File

@@ -1,21 +1,12 @@
BASENAME=pointprocesses
PYFILES=$(wildcard *.py)
PYPDFFILES=$(PYFILES:.py=.pdf)
GPTFILES=$(wildcard *.gpt)
GPTTEXFILES=$(GPTFILES:.gpt=.tex)
all: pdf slides thumbs
include ../../chapter.mk
# script:
pdf : $(BASENAME)-chapter.pdf
$(BASENAME)-chapter.pdf : $(BASENAME)-chapter.tex $(BASENAME).tex ../../header.tex $(GPTTEXFILES) $(PYPDFFILES)
CHAPTER=$$(( $$(sed -n -e '/contentsline {chapter}/{s/.*numberline {\([0123456789]*\)}.*/\1/; p}' $(BASENAME).aux) - 1 )); \
PAGE=$$(sed -n -e '/contentsline {chapter}/{s/.*numberline {.*}.*}{\(.*\)}{chapter.*/\1/; p}' $(BASENAME).aux); \
sed -i -e "s/setcounter{page}{.*}/setcounter{page}{$$PAGE}/; s/setcounter{chapter}{.*}/setcounter{chapter}{$$CHAPTER}/" $(BASENAME)-chapter.tex
pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true
pdf : chapter
# slides:
@@ -31,31 +22,15 @@ $(BASENAME)-handout.pdf: $(BASENAME)-slides.tex $(GPTTEXFILES)
pdfnup --nup 2x4 --no-landscape --paper a4paper --trim "-1cm -1cm -1cm -1cm" --outfile $@ thumbsfoils.pdf # 1-19
rm thumbsfoils.*
watchpdf :
while true; do ! make -q pdf && make pdf; sleep 0.5; done
watchslides :
while true; do ! make -q slides && make slides; sleep 0.5; done
# python plots:
$(PYPDFFILES) : %.pdf: %.py
python $<
# gnuplot plots:
$(GPTTEXFILES) : %.tex: %.gpt whitestyles.gp
gnuplot whitestyles.gp $<
epstopdf $*.eps
clean :
rm -f *~
rm -f $(BASENAME).aux $(BASENAME).log
rm -f $(BASENAME)-chapter.aux $(BASENAME)-chapter.log $(BASENAME)-chapter.out
clean : cleanchapter
rm -f $(BASENAME)-slides.aux $(BASENAME)-slides.log $(BASENAME)-slides.out $(BASENAME)-slides.toc $(BASENAME)-slides.nav $(BASENAME)-slides.snm $(BASENAME)-slides.vrb
rm -f $(PYPDFFILES) $(GPTTEXFILES)
cleanall : clean
rm -f $(BASENAME)-chapter.pdf $(BASENAME)-slides.pdf $(BASENAME)-handout.pdf
cleanall : clean cleanchapter
$(BASENAME)-slides.pdf $(BASENAME)-handout.pdf
help :

View File

@@ -98,4 +98,4 @@ plotisih(ax, isis(inhspikes))
plt.tight_layout()
plt.savefig('isihexamples.pdf')
plt.show()
plt.close()

View File

@@ -5,8 +5,10 @@
\lstset{inputpath=../code}
\graphicspath{{figures/}}
\setcounter{page}{111}
\setcounter{chapter}{7}
\typein[\pagenumber]{Number of first page}
\typein[\chapternumber]{Chapter number}
\setcounter{page}{\pagenumber}
\setcounter{chapter}{\chapternumber}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -14,5 +16,11 @@
\input{pointprocesses}
\section{TODO}
\begin{itemize}
\item Add spikeraster function
\item Multitrial firing rates
\end{itemize}
\end{document}

View File

@@ -2,7 +2,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Analyse von Spiketrains}
\determ{Aktionspotentiale} (\enterm{Spikes}) sind die Tr\"ager der
\determ[Aktionspotential]{Aktionspotentiale} (\enterm{spikes}) sind die Tr\"ager der
Information in Nervensystemen. Dabei ist in erster Linie nur der
Zeitpunkt des Auftretens eines Aktionspotentials von Bedeutung. Die
genaue Form des Aktionspotentials spielt keine oder nur eine
@@ -10,13 +10,13 @@ untergeordnete Rolle.
Nach etwas Vorverarbeitung haben elektrophysiologische Messungen
deshalb Listen von Spikezeitpunkten als Ergebniss --- sogenannte
\enterm{Spiketrains}. Diese Messungen k\"onnen wiederholt werden und
\enterm{spiketrains}. Diese Messungen k\"onnen wiederholt werden und
es ergeben sich mehrere \enterm{trials} von Spiketrains
(\figref{rasterexamplesfig}).
Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen
--- und deren Analyse f\"allt daher in das Gebiet der Statistik von
sogenannten \determ{Punktprozessen}.
sogenannten \determ[Punktprozess]{Punktprozessen}.
\begin{figure}[ht]
\includegraphics[width=1\textwidth]{rasterexamples}
@@ -25,7 +25,9 @@ sogenannten \determ{Punktprozessen}.
(homogener Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und
eines nicht-station\"aren Punktprozesses (perfect
integrate-and-fire Neuron getrieben mit Ohrnstein-Uhlenbeck
Rauschen mit Zeitkonstante $\tau=100$\,ms, rechts).}
Rauschen mit Zeitkonstante $\tau=100$\,ms, rechts). Jeder
vertikale Strich markiert den Zeitpunkt eines Ereignisses.
Jede Zeile zeigt die Ereignisse eines trials.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -79,11 +81,12 @@ Zeitpunkte der Ereignisse durch senkrechte Striche markiert werden.
Die Intervalle $T_i=t_{i+1}-t_i$ zwischen aufeinanderfolgenden
Ereignissen sind reelle, positive Zahlen. Bei Aktionspotentialen
heisen die Intervalle auch \enterm{Interspikeintervalle}. Deren Statistik
kann mit den \"ublichen Gr\"o{\ss}en beschrieben werden.
heisen die Intervalle auch \determ{Interspikeintervalle}
(\enterm{interspike intervals}). Deren Statistik kann mit den
\"ublichen Gr\"o{\ss}en beschrieben werden.
\begin{figure}[t]
\includegraphics[width=1\textwidth]{isihexamples}\hfill
\includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
\titlecaption{\label{isihexamplesfig}Interspikeintervall Histogramme}{der in
\figref{rasterexamplesfig} gezeigten Spikes.}
\end{figure}
@@ -104,21 +107,21 @@ kann mit den \"ublichen Gr\"o{\ss}en beschrieben werden.
\frac{1}{n}\sum\limits_{i=1}^n T_i$.
\item Standardabweichung der Intervalle: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
\rangle)^2 \rangle}$\vspace{1ex}
\item Variationskoeffizient (\enterm{coefficient of variation}): $CV_{ISI} =
\item \determ{Variationskoeffizient} (\enterm{coefficient of variation}): $CV_{ISI} =
\frac{\sigma_{ISI}}{\mu_{ISI}}$.
\item Diffusions Koeffizient: $D_{ISI} =
\item \determ{Diffusionskoeffizient} (\enterm{diffusion coefficient}): $D_{ISI} =
\frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
\end{itemize}
\begin{exercise}{isi_hist.m}{}
Schreibe eine Funktion \code{isi\_hist()}, die einen Vektor mit Interspikeintervallen
\begin{exercise}{isiHist.m}{}
Schreibe eine Funktion \code{isiHist()}, die einen Vektor mit Interspikeintervallen
entgegennimmt und daraus ein normiertes Histogramm der Interspikeintervalle
berechnet.
\end{exercise}
\begin{exercise}{plot_isi_hist.m}{}
\begin{exercise}{plotISIHist.m}{}
Schreibe eine Funktion, die die Histogrammdaten der Funktion
\code{isi\_hist()} entgegennimmt, um das Histogramm zu plotten. Im
\code{isiHist()} entgegennimmt, um das Histogramm zu plotten. Im
Plot sollen die Interspikeintervalle in Millisekunden aufgetragen
werden. Das Histogramm soll zus\"atzlich mit Mittelwert,
Standardabweichung und Variationskoeffizient der
@@ -139,9 +142,10 @@ sichtbar.
im Abstand des Lags $k$.}
\end{figure}
Solche Ab\"angigkeiten werden durch die serielle Korrelation der
Intervalle quantifiziert. Das ist der Korrelationskoeffizient
zwischen aufeinander folgenden Intervallen getrennt durch lag $k$:
Solche Ab\"angigkeiten werden durch die \determ{serielle
Korrelationen} (\enterm{serial correlations}) der Intervalle
quantifiziert. Das ist der \determ{Korrelationskoeffizient} zwischen
aufeinander folgenden Intervallen getrennt durch 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) \]
\"Ublicherweise wird die Korrelation $\rho_k$ gegen den Lag $k$
@@ -151,6 +155,7 @@ Intervalls mit sich selber).
\begin{exercise}{isiserialcorr.m}{}
Schreibe eine Funktion \code{isiserialcorr()}, die einen Vektor mit Interspikeintervallen
entgegennimmt und daraus die seriellen Korrelationen berechnet und plottet.
\pagebreak[4]
\end{exercise}
@@ -170,10 +175,10 @@ durch folgende Sch\"atzer charakterisiert werden:
\item Histogramm der counts $n_i$.
\item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$.
\item Varianz der Anzahl: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
\item Fano Faktor (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_n^2}{\mu_n}$.
\item \determ{Fano Faktor} (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_n^2}{\mu_n}$.
\end{itemize}
Insbesondere ist die mittlere Rate der Ereignisse $r$ (Spikes pro
Zeit, \determ{Feuerrate}) gemessen in Hertz
Zeit, \determ{Feuerrate}) gemessen in Hertz \sindex[term]{Feuerrate!mittlere Rate}
\begin{equation}
\label{firingrate}
r = \frac{\langle n \rangle}{W} \; .
@@ -209,18 +214,18 @@ u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine
\"ahnliche Rolle spielt bei Punktprozessen der \determ{Poisson
Prozess}.
Beim homogenen Poisson Prozess treten Ereignisse mit einer festen Rate
$\lambda=\text{const.}$ auf und sind unabh\"angig von der Zeit $t$ und
unabh\"angig von den Zeitpunkten fr\"uherer Ereignisse
(\figref{hompoissonfig}). Die Wahrscheinlichkeit zu irgendeiner Zeit
ein Ereigniss in einem kleinen Zeitfenster der Breite $\Delta t$ zu
bekommen ist
Beim \determ[Poisson Prozess!homogener]{homogenen Poisson Prozess}
treten Ereignisse mit einer festen Rate $\lambda=\text{const.}$ auf
und sind unabh\"angig von der Zeit $t$ und unabh\"angig von den
Zeitpunkten fr\"uherer Ereignisse (\figref{hompoissonfig}). Die
Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen
Zeitfenster der Breite $\Delta t$ zu bekommen ist
\begin{equation}
\label{hompoissonprob}
P = \lambda \cdot \Delta t \; .
\end{equation}
Beim inhomogenen Poisson Prozess h\"angt die Rate $\lambda$ von der
Zeit ab: $\lambda = \lambda(t)$.
Beim \determ[Poisson Prozess!inhomogener]{inhomogenen Poisson Prozess}
h\"angt die Rate $\lambda$ von der Zeit ab: $\lambda = \lambda(t)$.
\begin{exercise}{poissonspikes.m}{}
Schreibe eine Funktion \code{poissonspikes()}, die die Spikezeiten
@@ -253,14 +258,15 @@ Der homogene Poissonprozess hat folgende Eigenschaften:
\item Das mittlere Intervall ist $\mu_{ISI} = \frac{1}{\lambda}$ .
\item Die Varianz der Intervalle ist $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ .
\item Der Variationskoeffizient ist also immer $CV_{ISI} = 1$ .
\item Die seriellen Korrelationen $\rho_k =0$ f\"ur $k>0$, da das
Auftreten der Ereignisse unabh\"angig von der Vorgeschichte ist. Ein
solcher Prozess wird auch \determ{Erneuerungsprozess} genannt (\enterm{renewal
process}).
\item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der L\"ange W ist Poissonverteilt:
\item Die \determ[serielle Korrelationen]{seriellen Korrelationen}
$\rho_k =0$ f\"ur $k>0$, da das Auftreten der Ereignisse
unabh\"angig von der Vorgeschichte ist. Ein solcher Prozess wird
auch \determ{Erneuerungsprozess} genannt (\enterm{renewal process}).
\item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der
L\"ange W ist \determ[Poisson-Verteilung]{Poissonverteilt}:
\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]
(\figref{hompoissoncountfig})
\item Der Fano Faktor ist immer $F=1$ .
\item Der \determ{Fano Faktor} ist immer $F=1$ .
\end{itemize}
\begin{exercise}{hompoissonspikes.m}{}
@@ -299,13 +305,11 @@ Abbildung \ref{psthfig} n\"aher erl\"autert.
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{firingrates}
\titlecaption{Verschiedene Methoden die zeitabh\"angige Feuerrate
zu bestimmen.}{\textbf{A)} Rasterplot einer einzelnen neuronalen
Antwort. Jeder vertikale Strich notiert den Zeitpunkt eines
Aktionspotentials. \textbf{B)} Feurerrate aus der instantanen
Feuerrate bestimmt. \textbf{C)} klassisches PSTH mit der Binning
Methode. \textbf{D)} Feuerrate durch Faltung mit einem Gauss Kern
bestimmt.}\label{psthfig}
\titlecaption{Bestimmung der zeitabh\"angigen
Feuerrate.}{\textbf{A)} Rasterplot eines Spiketrains. \textbf{B)}
Feurerrate aus der instantanen Feuerrate bestimmt. \textbf{C)}
klassisches PSTH mit der Binning Methode. \textbf{D)} Feuerrate
durch Faltung mit einem Gauss Kern bestimmt.}\label{psthfig}
\end{figure}
@@ -314,26 +318,27 @@ Abbildung \ref{psthfig} n\"aher erl\"autert.
\begin{figure}[tp]
\includegraphics[width=\columnwidth]{isimethod}
\titlecaption{Instantane Feuerrate.}{Skizze eines Spiketrains
(oben). Jeder vertikale Strich notiert den Zeitpunkt eines
Aktionspotentials. Die Pfeile zwischen aufeinanderfolgenden
(oben). Die Pfeile zwischen aufeinanderfolgenden
Aktionspotentialen mit den Zahlen in Millisekunden illustrieren
die Interspikeintervalle. Der Kehrwert des Interspikeintervalle
ergibt die instantane Feuerrate.}\label{instrate}
\end{figure}
Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist
die sogenannte \determ{instantane Feuerrate}. Dabei wird die Feuerrate
aus dem Kehrwert der Interspikeintervalle, der Zeit zwischen zwei
aufeinander folgenden Aktionspotentialen (\figref{instrate} A),
bestimmt. Die abgesch\"atzte Feuerrate (\figref{instrate} B) ist
g\"ultig f\"ur das gesammte Interspikeintervall. Diese Methode hat den
Vorteil, dass sie sehr einfach zu berechnen ist und keine Annahme
\"uber eine relevante Zeitskala (der Kodierung oder des
Auslesemechanismus der postsynaptischen Zelle) macht. $r(t)$ ist
allerdings keine kontinuierliche Funktion, die Spr\"unge in der
Feuerrate k\"onnen f\"ur manche Analysen nachteilig sein. Au{\ss}erdem
wird die Feuerrate nie gleich Null, auch wenn lange keine Aktionspotentiale
generiert wurden.
die sogenannte \determ[Feuerrate!instantane]{instantane Feuerrate}
(\enterm[firing rate!instantaneous]{instantaneous firing rate}). Dabei
wird die Feuerrate aus dem Kehrwert der Interspikeintervalle, der Zeit
zwischen zwei aufeinander folgenden Aktionspotentialen
(\figref{instrate} A), bestimmt. Die abgesch\"atzte Feuerrate
(\figref{instrate} B) ist g\"ultig f\"ur das gesammte
Interspikeintervall. Diese Methode hat den Vorteil, dass sie sehr
einfach zu berechnen ist und keine Annahme \"uber eine relevante
Zeitskala (der Kodierung oder des Auslesemechanismus der
postsynaptischen Zelle) macht. $r(t)$ ist allerdings keine
kontinuierliche Funktion, die Spr\"unge in der Feuerrate k\"onnen
f\"ur manche Analysen nachteilig sein. Au{\ss}erdem wird die Feuerrate
nie gleich Null, auch wenn lange keine Aktionspotentiale generiert
wurden.
\begin{exercise}{instantaneousRate.m}{}
Implementiere die Absch\"atzung der Feuerrate auf Basis der
@@ -345,9 +350,10 @@ generiert wurden.
W\"ahrend die Instantane Rate den Kehrwert der Zeit von einem bis zum
n\"achsten Aktionspotential misst, sch\"atzt das sogenannte
\determ{Peri-Stimulus-Zeit-Histogramm} (\enterm{peri stimulus time
histogram}, PSTH) die Wahrscheinlichkeit ab, zu einem Zeitpunkt
Aktionspotentiale anzutreffen. Es wird versucht die mittlere Rate \eqnref{firingrate}
im Grenzwert kleiner Beobachtungszeiten abzusch\"atzen:
histogram}, \determ[PSTH|see{Peri-Stimulus-Zeit-Histogramm}]{PSTH})
die Wahrscheinlichkeit ab, zu einem Zeitpunkt Aktionspotentiale
anzutreffen. Es wird versucht die mittlere Rate \eqnref{firingrate} im
Grenzwert kleiner Beobachtungszeiten abzusch\"atzen:
\begin{equation}
\label{psthrate}
r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; ,
@@ -377,9 +383,9 @@ Bei der Binning-Methode wird die Zeitachse in gleichm\"aßige
Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in
die jeweiligen Bins fallen, gez\"ahlt (\figref{binpsth} A). Um diese
Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der Binweite
normiert werden. Das ist fast so, wie beim Absch\"atzen einer
Wahrscheinlichkeitsdichte. Es kann auch die \code{hist} Funktion zur
Bestimmung des PSTHs verwendet werden.
normiert werden. Das ist \"aquivalent zur Absch\"atzung einer
Wahrscheinlichkeitsdichte. Es kann auch die \code{hist()} Funktion zur
Bestimmung des PSTHs verwendet werden. \sindex[term]{Feuerrate!Binningmethode}
Die bestimmte Feuerrate gilt f\"ur das gesamte Bin (\figref{binpsth}
B). Das so berechnete PSTH hat wiederum eine stufige Form, die von der
@@ -390,6 +396,7 @@ vorkommen k\"onnen nicht aufgl\"ost werden. Mit der Wahl der Binweite
wird somit eine Annahme \"uber die relevante Zeitskala des Spiketrains
gemacht.
\pagebreak[4]
\begin{exercise}{binnedRate.m}{}
Implementiere die Absch\"atzung der Feuerrate mit der ``binning''
Methode. Plotte das PSTH.
@@ -421,7 +428,7 @@ wobei $\omega(\tau)$ der Filterkern und $\rho(t)$ die bin\"are Antwort
ist. Bildlich geprochen wird jede 1 in $\rho(t)$ durch den Filterkern
ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert
wurde (Integral gleich Eins), ergibt sich die Feuerrate direkt aus der
\"Uberlagerung der Kerne (Abb. \ref{convrate} B).
\"Uberlagerung der Kerne (Abb. \ref{convrate} B). \sindex[term]{Feuerrate!Faltungsmethode}
Die Faltungsmethode f\"uhrt, anders als die anderen Methoden, zu einer
stetigen Funktion was insbesondere f\"ur spektrale Analysen von
@@ -430,6 +437,7 @@ Binweite, die zeitliche Aufl\"osung von $r(t)$. Die Breite des Kerns
macht also auch wieder eine Annahme \"uber die relevante Zeitskala des
Spiketrains.
\pagebreak[4]
\begin{exercise}{convolutionRate.m}{}
Verwende die Faltungsmethode um die Feuerrate zu bestimmen. Plotte
das Ergebnis.
@@ -438,8 +446,9 @@ Spiketrains.
\section{Spike-triggered Average}
Die graphischer Darstellung der Feuerrate allein reicht nicht aus um
den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu
analysieren. Eine Methode um mehr \"uber diesen Zusammenhang zu erfahren,
ist der \enterm{Spike-triggered average} (STA). Der STA
analysieren. Eine Methode um mehr \"uber diesen Zusammenhang zu
erfahren, ist der \enterm{spike-triggered average}
(\enterm[STA|see{spike-triggered average}]{STA}). Der STA
\begin{equation}
STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau)
\end{equation}
@@ -477,17 +486,20 @@ antworten.
Der STA kann auch dazu benutzt werden, aus den Antworten der Zelle den
Stimulus zu rekonstruieren (\figref{stafig} B). Bei der
\determ{invertierten Rekonstruktion} wird die Zellantwort mit dem STA
verfaltet.
\determ[invertierte Rekonstruktion]{invertierten Rekonstruktion} wird
die Zellantwort mit dem STA verfaltet.
\begin{exercise}{spikeTriggeredAverage.m}{}
Implementiere eine Funktion, die den STA ermittelt. Verwende dazu
den Datensatz \codeterm{sta\_data.mat}. Die Funktion sollte folgende
den Datensatz \file{sta\_data.mat}. Die Funktion sollte folgende
R\"uckgabewerte haben:
\vspace{-1ex}
\begin{itemize}
\setlength{\itemsep}{0ex}
\item den Spike-Triggered-Average.
\item die Standardabweichung der individuellen STAs.
\item die Anzahl Aktionspotentiale, die dem STA zugrunde liegen.
\item die Anzahl Aktionspotentiale, die zur Berechnung des STA verwendet wurden.
\vspace{-2ex}
\end{itemize}
\end{exercise}
@@ -495,7 +507,5 @@ verfaltet.
Rekonstruiere den Stimulus mithilfe des STA und der Spike
Zeiten. Die Funktion soll Vektor als R\"uckgabewert haben, der
genauso gro{\ss} ist wie der Originalstimulus aus der Datei
\codeterm{sta\_data.mat}.
\file{sta\_data.mat}.
\end{exercise}

View File

@@ -1,4 +1,4 @@
set term epslatex size 11.4cm, 7cm
set term epslatex size 11.4cm, 6.5cm
set out 'pointprocessscetch.tex'
set border 0

View File

@@ -1,7 +1,7 @@
%!PS-Adobe-2.0 EPSF-2.0
%%Title: pointprocessscetchA.tex
%%Creator: gnuplot 4.6 patchlevel 4
%%CreationDate: Tue Nov 3 17:29:16 2015
%%CreationDate: Sat Nov 28 12:01:31 2015
%%DocumentFonts:
%%BoundingBox: 50 50 373 135
%%EndComments
@@ -430,10 +430,10 @@ SDict begin [
/Title (pointprocessscetchA.tex)
/Subject (gnuplot plot)
/Creator (gnuplot 4.6 patchlevel 4)
/Author (grewe)
/Author (jan)
% /Producer (gnuplot)
% /Keywords ()
/CreationDate (Tue Nov 3 17:29:16 2015)
/CreationDate (Sat Nov 28 12:01:31 2015)
/DOCINFO pdfmark
end
} ifelse

View File

@@ -1,7 +1,7 @@
%!PS-Adobe-2.0 EPSF-2.0
%%Title: pointprocessscetchB.tex
%%Creator: gnuplot 4.6 patchlevel 4
%%CreationDate: Tue Nov 3 18:24:39 2015
%%CreationDate: Sat Nov 28 12:01:32 2015
%%DocumentFonts:
%%BoundingBox: 50 50 373 237
%%EndComments
@@ -430,10 +430,10 @@ SDict begin [
/Title (pointprocessscetchB.tex)
/Subject (gnuplot plot)
/Creator (gnuplot 4.6 patchlevel 4)
/Author (grewe)
/Author (jan)
% /Producer (gnuplot)
% /Keywords ()
/CreationDate (Tue Nov 3 18:24:39 2015)
/CreationDate (Sat Nov 28 12:01:32 2015)
/DOCINFO pdfmark
end
} ifelse

View File

@@ -83,4 +83,4 @@ ax.eventplot(inhspikes, colors=[[0, 0, 0]], linelength=0.8)
plt.tight_layout()
plt.savefig('rasterexamples.pdf')
plt.show()
plt.close()

View File

@@ -38,18 +38,22 @@ def reconstruct_stimulus(spike_times, sta, stimulus, t_max=30., dt=1e-4):
def plot_results(sta_time, st_average, stim_time, s_est, stimulus, duration, dt):
plt.xkcd()
sta_ax = plt.subplot2grid((1, 3), (0, 0), rowspan=1, colspan=1)
stim_ax = plt.subplot2grid((1, 3), (0, 1), rowspan=1, colspan=2)
fig = plt.gcf()
fig.set_size_inches(15, 5)
fig.subplots_adjust(left=0.075, bottom=0.12, top=0.92, right=0.975)
fig.set_size_inches(8, 3)
fig.subplots_adjust(left=0.08, bottom=0.15, top=0.9, right=0.975)
fig.set_facecolor("white")
sta_ax.plot(sta_time * 1000, st_average, color="dodgerblue", lw=2.)
sta_ax.set_xlabel("time [ms]", fontsize=12)
sta_ax.set_ylabel("stimulus", fontsize=12)
sta_ax.set_xlim([-50, 50])
sta_ax.plot(sta_time * 1000, st_average, color="#FF9900", lw=2.)
sta_ax.set_xlabel("Time (ms)")
sta_ax.set_ylabel("Stimulus")
sta_ax.set_xlim(-40, 20)
sta_ax.set_xticks(np.arange(-40, 21, 20))
sta_ax.set_ylim(-0.1, 0.2)
sta_ax.set_yticks(np.arange(-0.1, 0.21, 0.1))
# sta_ax.xaxis.grid('off')
sta_ax.spines["right"].set_visible(False)
sta_ax.spines["top"].set_visible(False)
@@ -58,22 +62,31 @@ def plot_results(sta_time, st_average, stim_time, s_est, stimulus, duration, dt)
sta_ax.spines["bottom"].set_linewidth(2.0)
sta_ax.spines["left"].set_linewidth(2.0)
sta_ax.tick_params(direction="out", width=2.0)
ylim = sta_ax.get_ylim()
xlim = sta_ax.get_xlim()
sta_ax.plot(list(xlim), [0., 0.], zorder=1, color='darkgray', ls='--', lw=0.75)
sta_ax.plot([0., 0.], list(ylim), zorder=1, color='darkgray', ls='--', lw=0.75)
sta_ax.plot(list(xlim), [0., 0.], zorder=1, color='darkgray', ls='--', lw=1)
sta_ax.plot([0., 0.], list(ylim), zorder=1, color='darkgray', ls='--', lw=1)
sta_ax.set_xlim(list(xlim))
sta_ax.set_ylim(list(ylim))
sta_ax.text(-0.225, 1.05, "A", transform=sta_ax.transAxes, size=14)
sta_ax.annotate('Time of\nspike',
xy=(0, 0.18), xycoords='data',
xytext=(-35, 0.19), textcoords='data', ha='left',
arrowprops=dict(arrowstyle="->", relpos=(1.0,0.5),
connectionstyle="angle3,angleA=0,angleB=-70") )
sta_ax.annotate('STA',
xy=(-10, 0.05), xycoords='data',
xytext=(-33, 0.09), textcoords='data', ha='left',
arrowprops=dict(arrowstyle="->", relpos=(1.0,0.0),
connectionstyle="angle3,angleA=60,angleB=-40") )
#sta_ax.text(-0.25, 1.04, "A", transform=sta_ax.transAxes, size=24)
stim_ax.plot(stim_time * 1000, stimulus[:,1], label='stimulus', color='dodgerblue', lw=2.)
stim_ax.plot(stim_time * 1000, s_est, label='reconstruction', color='red', lw=2)
stim_ax.set_xlabel('time[ms]', fontsize=12)
stim_ax.set_xlim([0.0, 250])
stim_ax.plot(stim_time * 1000, stimulus[:,1], label='stimulus', color='#0000FF', lw=2.)
stim_ax.plot(stim_time * 1000, s_est, label='reconstruction', color='#FF9900', lw=2)
stim_ax.set_xlabel('Time (ms)')
stim_ax.set_xlim(0.0, 200)
stim_ax.set_ylim([-1., 1.])
stim_ax.legend()
stim_ax.plot([0.0, 250], [0., 0.], color="darkgray", lw=0.75, ls='--', zorder=1)
stim_ax.legend(loc=(0.3, 0.85), frameon=False, fontsize=12)
stim_ax.plot([0.0, 250], [0., 0.], color="darkgray", lw=1, ls='--', zorder=1)
stim_ax.spines["right"].set_visible(False)
stim_ax.spines["top"].set_visible(False)
stim_ax.yaxis.set_ticks_position('left')
@@ -81,8 +94,9 @@ def plot_results(sta_time, st_average, stim_time, s_est, stimulus, duration, dt)
stim_ax.spines["bottom"].set_linewidth(2.0)
stim_ax.spines["left"].set_linewidth(2.0)
stim_ax.tick_params(direction="out", width=2.0)
stim_ax.text(-0.075, 1.05, "B", transform=stim_ax.transAxes, size=14)
#stim_ax.text(-0.1, 1.04, "B", transform=stim_ax.transAxes, size=24)
fig.tight_layout()
fig.savefig("sta.pdf")
plt.close()