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

Merge branch 'master' of https://whale.am28.uni-tuebingen.de/git/teaching/scientificComputing

Browse Source
This commit is contained in:
Jan Grewe 2019-01-13 11:24:20 +01:00
commit 263c08e09c
65 changed files with 1449 additions and 396 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

192
likelihood/exercises/exercises01-de.tex Normal file
View File

@ -0,0 +1,192 @@
\documentclass[12pt,a4paper,pdftex]{exam}
\usepackage[german]{babel}
\usepackage{pslatex}
\usepackage[mediumspace,mediumqspace,Gray]{SIunits} % \ohm, \micro
\usepackage{xcolor}
\usepackage{graphicx}
\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue]{hyperref}
%%%%% layout %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{headandfoot}
\ifprintanswers
\newcommand{\stitle}{: L\"osungen}
\else
\newcommand{\stitle}{}
\fi
\header{{\bfseries\large \"Ubung 4\stitle}}{{\bfseries\large Statistik}}{{\bfseries\large 26. Oktober, 2015}}
\firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email:
jan.benda@uni-tuebingen.de}
\runningfooter{}{\thepage}{}
\setlength{\baselineskip}{15pt}
\setlength{\parindent}{0.0cm}
\setlength{\parskip}{0.3cm}
\renewcommand{\baselinestretch}{1.15}
%%%%% listings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{listings}
\lstset{
language=Matlab,
basicstyle=\ttfamily\footnotesize,
numbers=left,
numberstyle=\tiny,
title=\lstname,
showstringspaces=false,
commentstyle=\itshape\color{darkgray},
breaklines=true,
breakautoindent=true,
columns=flexible,
frame=single,
xleftmargin=1em,
xrightmargin=1em,
aboveskip=10pt
}
%%%%% math stuff: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{bm}
\usepackage{dsfont}
\newcommand{\naZ}{\mathds{N}}
\newcommand{\gaZ}{\mathds{Z}}
\newcommand{\raZ}{\mathds{Q}}
\newcommand{\reZ}{\mathds{R}}
\newcommand{\reZp}{\mathds{R^+}}
\newcommand{\reZpN}{\mathds{R^+_0}}
\newcommand{\koZ}{\mathds{C}}
%%%%% page breaks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\continue}{\ifprintanswers%
\else
\vfill\hspace*{\fill}$\rightarrow$\newpage%
\fi}
\newcommand{\continuepage}{\ifprintanswers%
\newpage
\else
\vfill\hspace*{\fill}$\rightarrow$\newpage%
\fi}
\newcommand{\newsolutionpage}{\ifprintanswers%
\newpage%
\else
\fi}
%%%%% new commands %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\qt}[1]{\textbf{#1}\\}
\newcommand{\pref}[1]{(\ref{#1})}
\newcommand{\extra}{--- Zusatzaufgabe ---\ \mbox{}}
\newcommand{\code}[1]{\texttt{#1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\input{instructions}
\begin{questions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Maximum Likelihood der Standardabweichung}
Wir wollen uns die Likelihood und die Log-Likelihood am Beispiel der
Absch\"atzung der Standardabweichung verdeutlichen.
\begin{parts}
\part Ziehe $n=50$ normalverteilte Zufallsvariablen mit Mittelwert $\mu=3$
und einer Standardabweichung $\sigma=2$.
\part
Plotte die Likelihood (aus dem Produkt der Wahrscheinlichkeiten) und
die Log-Likelihood (aus der Summe der logarithmierten
Wahrscheinlichkeiten) f\"ur die Standardabweichung als Parameter. Vergleiche die
Position der Maxima mit der aus den Daten berechneten Standardabweichung.
\part
Erh\"ohe $n$ auf 1000. Was passiert mit der Likelihood, was mit der Log-Likelihood? Warum?
\end{parts}
\begin{solution}
\lstinputlisting{mlestd.m}
\includegraphics[width=1\textwidth]{mlestd}
\end{solution}
\continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Maximum-Likelihood-Sch\"atzer einer Ursprungsgeraden}
In der Vorlesung haben wir folgende Formel f\"ur die Maximum-Likelihood
Absch\"atzung der Steigung $\theta$ einer Ursprungsgeraden durch $n$ Datenpunkte $(x_i|y_i)$ mit Standardabweichung $\sigma_i$ hergeleitet:
\[\theta = \frac{\sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2}}{ \sum_{i=1}^n
\frac{x_i^2}{\sigma_i^2}} \]
\begin{parts}
\part \label{mleslopefunc} Schreibe eine Funktion, die in einem $x$ und einem
$y$ Vektor die Datenpaare \"uberreicht bekommt und die Steigung der
Ursprungsgeraden, die die Likelihood maximiert, zur\"uckgibt
($\sigma=\text{const}$).
\part
Schreibe ein Skript, das Datenpaare erzeugt, die um eine
Ursprungsgerade mit vorgegebener Steigung streuen. Berechne mit der
Funktion aus \pref{mleslopefunc} die Steigung aus den Daten,
vergleiche mit der wahren Steigung, und plotte die urspr\"ungliche
sowie die gefittete Gerade zusammen mit den Daten.
\part
Ver\"andere die Anzahl der Datenpunkte, die Steigung, sowie die
Streuung der Daten um die Gerade.
\end{parts}
\begin{solution}
\lstinputlisting{mleslope.m}
\lstinputlisting{mlepropfit.m}
\includegraphics[width=1\textwidth]{mlepropfit}
\end{solution}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Maximum-Likelihood-Sch\"atzer einer Wahrscheinlichkeitsdichtefunktion}
Verschiedene Wahrscheinlichkeitsdichtefunktionen haben Parameter, die
nicht so einfach wie der Mittelwert und die Standardabweichung einer
Normalverteilung direkt aus den Daten berechnet werden k\"onnen. Solche Parameter
m\"ussen dann aus den Daten mit der Maximum-Likelihood-Methode gefittet werden.
Um dies zu veranschaulichen ziehen wir uns diesmal nicht normalverteilte Zufallszahlen, sondern Zufallszahlen aus der Gamma-Verteilung.
\begin{parts}
\part
Finde heraus welche \code{matlab} Funktion die
Wahrscheinlichkeitsdichtefunktion (probability density function) der
Gamma-Verteilung berechnet.
\part
Plotte mit Hilfe dieser Funktion die Wahrscheinlichkeitsdichtefunktion
der Gamma-Verteilung f\"ur verschiedene Werte des (positiven) ``shape'' Parameters.
Den ``scale'' Parameter setzen wir auf Eins.
\part
Finde heraus mit welcher Funktion Gammaverteilte Zufallszahlen in
\code{matlab} gezogen werden k\"onnen. Erzeuge mit dieser Funktion
50 Zufallszahlen mit einem der oben geplotteten ``shape'' Parameter.
\part
Berechne und plotte ein normiertes Histogramm dieser Zufallszahlen.
\part
Finde heraus mit welcher \code{matlab}-Funktion eine beliebige
Verteilung (``distribution'') an die Zufallszahlen nach der
Maximum-Likelihood Methode gefittet werden kann. Wie wird diese
Funktion benutzt, um die Gammaverteilung an die Daten zu fitten?
\part
Bestimme mit dieser Funktion die Parameter der Gammaverteilung aus
den Zufallszahlen.
\part
Plotte anschlie{\ss}end die Gammaverteilung mit den gefitteten
Parametern.
\end{parts}
\begin{solution}
\lstinputlisting{mlepdffit.m}
\includegraphics[width=1\textwidth]{mlepdffit}
\end{solution}
\end{questions}
\end{document}

175
likelihood/exercises/exercises01.tex
View File

@ -11,11 +11,11 @@
\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{headandfoot}
\ifprintanswers
\newcommand{\stitle}{: L\"osungen}
\newcommand{\stitle}{: Solutions}
\else
\newcommand{\stitle}{}
\fi
\header{{\bfseries\large \"Ubung 4\stitle}}{{\bfseries\large Statistik}}{{\bfseries\large 26. Oktober, 2015}}
\header{{\bfseries\large Exercise 12\stitle}}{{\bfseries\large Maximum likelihood}}{{\bfseries\large January 7th, 2019}}
\firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email:
jan.benda@uni-tuebingen.de}
\runningfooter{}{\thepage}{}
@ -89,98 +89,119 @@ jan.benda@uni-tuebingen.de}
\begin{questions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Maximum Likelihood der Standardabweichung}
Wir wollen uns die Likelihood und die Log-Likelihood am Beispiel der
Absch\"atzung der Standardabweichung verdeutlichen.
\question \qt{Maximum likelihood of the standard deviation}
Let's compute the likelihood and the log-likelihood for the estimation
of the standard deviation.
\begin{parts}
\part Ziehe $n=50$ normalverteilte Zufallsvariablen mit Mittelwert $\mu=3$
und einer Standardabweichung $\sigma=2$.
\part Draw $n=50$ random numbers from a normal distribution with
mean $\mu=3$ and standard deviation $\sigma=2$.
\part
Plotte die Likelihood (aus dem Produkt der Wahrscheinlichkeiten) und
die Log-Likelihood (aus der Summe der logarithmierten
Wahrscheinlichkeiten) f\"ur die Standardabweichung als Parameter. Vergleiche die
Position der Maxima mit der aus den Daten berechneten Standardabweichung.
\part Plot the likelihood (computed as the product of probabilities)
and the log-likelihood (sum of the logarithms of the probabilities)
as a function of the standard deviation. Compare the position of the
maxima with the standard deviation that you compute directly from
the data.
\part
Erh\"ohe $n$ auf 1000. Was passiert mit der Likelihood, was mit der Log-Likelihood? Warum?
\part Increase $n$ to 1000. What happens to the likelihood, what
happens to the log-likelihood? Why?
\end{parts}
\begin{solution}
\lstinputlisting{mlestd.m}
\includegraphics[width=1\textwidth]{mlestd}
\includegraphics[width=1\textwidth]{mlestd}\\
The more data the smaller the product of the probabilities ($\approx
p^n$ with $0 \le p < 1$) and the smaller the sum of the logarithms
of the probabilities ($\approx n\log p$, note that $\log p < 0$).
The product eventually gets smaller than the precision of the
floating point numbers support. Therefore for $n=1000$ the products
becomes zero. Using the logarithm avoids this numerical problem.
\end{solution}
\continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Maximum-Likelihood-Sch\"atzer einer Ursprungsgeraden}
In der Vorlesung haben wir folgende Formel f\"ur die Maximum-Likelihood
Absch\"atzung der Steigung $\theta$ einer Ursprungsgeraden durch $n$ Datenpunkte $(x_i|y_i)$ mit Standardabweichung $\sigma_i$ hergeleitet:
\[\theta = \frac{\sum_{i=1}^n \frac{x_iy_i}{\sigma_i^2}}{ \sum_{i=1}^n
\question \qt{Maximum-likelihood estimator of a line through the origin}
In the lecture we derived the following equation for an
maximum-likelihood estimate of the slope $\theta$ of a straight line
through the origin fitted to $n$ pairs of data values $(x_i|y_i)$ with
standard deviation $\sigma_i$:
\[\theta = \frac{\sum_{i=1}^n \frac{x_i y_i}{\sigma_i^2}}{ \sum_{i=1}^n
\frac{x_i^2}{\sigma_i^2}} \]
\begin{parts}
\part \label{mleslopefunc} Schreibe eine Funktion, die in einem $x$ und einem
$y$ Vektor die Datenpaare \"uberreicht bekommt und die Steigung der
Ursprungsgeraden, die die Likelihood maximiert, zur\"uckgibt
($\sigma=\text{const}$).
\part
Schreibe ein Skript, das Datenpaare erzeugt, die um eine
Ursprungsgerade mit vorgegebener Steigung streuen. Berechne mit der
Funktion aus \pref{mleslopefunc} die Steigung aus den Daten,
vergleiche mit der wahren Steigung, und plotte die urspr\"ungliche
sowie die gefittete Gerade zusammen mit den Daten.
\part
Ver\"andere die Anzahl der Datenpunkte, die Steigung, sowie die
Streuung der Daten um die Gerade.
\part \label{mleslopefunc} Write a function that takes two vectors
$x$ and $y$ containing the data pairs and returns the slope,
computed according to this equation. For simplicity we assume
$\sigma_i=\sigma$ for all $1 \le i \le n$. How does this simplify
the equation for the slope?
\begin{solution}
\lstinputlisting{mleslope.m}
\end{solution}
\part Write a script that generates data pairs that scatter around a
line through the origin with a given slope. Use the function from
\pref{mleslopefunc} to compute the slope from the generated data.
Compare the computed slope with the true slope that has been used to
generate the data. Plot the data togehther with the line from which
the data were generated and the maximum-likelihood fit.
\begin{solution}
\lstinputlisting{mlepropfit.m}
\includegraphics[width=1\textwidth]{mlepropfit}
\end{solution}
\part \label{mleslopecomp} Vary the number of data pairs, the slope,
as well as the variance of the data points around the true
line. Under which conditions is the maximum-likelihood estimation of
the slope closer to the true slope?
\part To answer \pref{mleslopecomp} more precisely, generate for
each condition let's say 1000 data sets and plot a histogram of the
estimated slopes. How does the histogram, its mean and standard
deviation relate to the true slope?
\end{parts}
\begin{solution}
\lstinputlisting{mleslope.m}
\lstinputlisting{mlepropfit.m}
\includegraphics[width=1\textwidth]{mlepropfit}
\lstinputlisting{mlepropest.m}
\includegraphics[width=1\textwidth]{mlepropest}\\
The estimated slopes are centered around the true slope. The
standard deviation of the estimated slopes gets smaller for larger
$n$ and less noise in the data.
\end{solution}
\continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\question \qt{Maximum-Likelihood-Sch\"atzer einer Wahrscheinlichkeitsdichtefunktion}
Verschiedene Wahrscheinlichkeitsdichtefunktionen haben Parameter, die
nicht so einfach wie der Mittelwert und die Standardabweichung einer
Normalverteilung direkt aus den Daten berechnet werden k\"onnen. Solche Parameter
m\"ussen dann aus den Daten mit der Maximum-Likelihood-Methode gefittet werden.
Um dies zu veranschaulichen ziehen wir uns diesmal nicht normalverteilte Zufallszahlen, sondern Zufallszahlen aus der Gamma-Verteilung.
\question \qt{Maximum-likelihood-estimation of a probability-density function}
Many probability-density functions have parameters that cannot be
computed directly from the data, like, for example, the mean of
normally-distributed data. Such parameter need to be estimated by
means of the maximum-likelihood from the data.
Let us demonstrate this approach by means of data that are drawn from a
gamma distribution,
\begin{parts}
\part
Finde heraus welche \code{matlab} Funktion die
Wahrscheinlichkeitsdichtefunktion (probability density function) der
Gamma-Verteilung berechnet.
\part
Plotte mit Hilfe dieser Funktion die Wahrscheinlichkeitsdichtefunktion
der Gamma-Verteilung f\"ur verschiedene Werte des (positiven) ``shape'' Parameters.
Den ``scale'' Parameter setzen wir auf Eins.
\part
Finde heraus mit welcher Funktion Gammaverteilte Zufallszahlen in
\code{matlab} gezogen werden k\"onnen. Erzeuge mit dieser Funktion
50 Zufallszahlen mit einem der oben geplotteten ``shape'' Parameter.
\part
Berechne und plotte ein normiertes Histogramm dieser Zufallszahlen.
\part
Finde heraus mit welcher \code{matlab}-Funktion eine beliebige
Verteilung (``distribution'') an die Zufallszahlen nach der
Maximum-Likelihood Methode gefittet werden kann. Wie wird diese
Funktion benutzt, um die Gammaverteilung an die Daten zu fitten?
\part
Bestimme mit dieser Funktion die Parameter der Gammaverteilung aus
den Zufallszahlen.
\part
Plotte anschlie{\ss}end die Gammaverteilung mit den gefitteten
Parametern.
\part Find out which \code{matlab} function computes the
probability-density function of the gamma distribution.
\part \label{gammaplot} Use this function to plot the
probability-density function of the gamma distribution for various
values of the (positive) ``shape'' parameter. Wet set the ``scale''
parameter to one.
\part Find out which \code{matlab} function generates random numbers
that are distributed according to a gamma distribution. Generate
with this function 50 random numbers using one of the values of the
``shape'' parameter used in \pref{gammaplot}.
\part Compute and plot a properly normalized histogram of these
random numbers.
\part Find out which \code{matlab} function fit a distribution to a
vector of random numbers according to the maximum-likelihood method.
How do you need to use this function in order to fit a gamma
distribution to the data?
\part Estimate with this function the parameter of the gamma
distribution used to generate the data.
\part Finally, plot the fitted gamma distribution on top of the
normalized histogram of the data.
\end{parts}
\begin{solution}
\lstinputlisting{mlepdffit.m}

68
likelihood/exercises/instructions.tex
View File

@ -1,41 +1,41 @@
\vspace*{-6.5ex}
\vspace*{-8ex}
\begin{center}
\textbf{\Large Einf\"uhrung in die wissenschaftliche Datenverarbeitung}\\[1ex]
\textbf{\Large Introduction to scientific computing}\\[1ex]
{\large Jan Grewe, Jan Benda}\\[-3ex]
Abteilung Neuroethologie \hfill --- \hfill Institut f\"ur Neurobiologie \hfill --- \hfill \includegraphics[width=0.28\textwidth]{UT_WBMW_Black_RGB} \\
Neuroethology lab \hfill --- \hfill Institute for Neurobiology \hfill --- \hfill \includegraphics[width=0.28\textwidth]{UT_WBMW_Black_RGB} \\
\end{center}
\ifprintanswers%
\else
% \ifprintanswers%
% \else
% Die folgenden Aufgaben dienen der Wiederholung, \"Ubung und
% Selbstkontrolle und sollten eigenst\"andig bearbeitet und gel\"ost
% werden. Die L\"osung soll in Form eines einzelnen Skriptes (m-files)
% im ILIAS hochgeladen werden. Jede Aufgabe sollte in einer eigenen
% ``Zelle'' gel\"ost sein. Die Zellen \textbf{m\"ussen} unabh\"angig
% voneinander ausf\"uhrbar sein. Das Skript sollte nach dem Muster:
% ``variablen\_datentypen\_\{nachname\}.m'' benannt werden
% (z.B. variablen\_datentypen\_mueller.m).
% % Die folgenden Aufgaben dienen der Wiederholung, \"Ubung und
% % Selbstkontrolle und sollten eigenst\"andig bearbeitet und gel\"ost
% % werden. Die L\"osung soll in Form eines einzelnen Skriptes (m-files)
% % im ILIAS hochgeladen werden. Jede Aufgabe sollte in einer eigenen
% % ``Zelle'' gel\"ost sein. Die Zellen \textbf{m\"ussen} unabh\"angig
% % voneinander ausf\"uhrbar sein. Das Skript sollte nach dem Muster:
% % ``variablen\_datentypen\_\{nachname\}.m'' benannt werden
% % (z.B. variablen\_datentypen\_mueller.m).
\begin{itemize}
\item \"Uberzeuge dich von jeder einzelnen Zeile deines Codes, dass
sie auch wirklich das macht, was sie machen soll! Teste dies mit
kleinen Beispielen direkt in der Kommandozeile.
\item Versuche die L\"osungen der Aufgaben m\"oglichst in
sinnvolle kleine Funktionen herunterzubrechen.
Sobald etwas \"ahnliches mehr als einmal berechnet werden soll,
lohnt es sich eine Funktion daraus zu schreiben!
\item Teste rechenintensive \code{for} Schleifen, Vektoren, Matrizen
zuerst mit einer kleinen Anzahl von Wiederholungen oder kleiner
Gr\"o{\ss}e, und benutze erst am Ende, wenn alles \"uberpr\"uft
ist, eine gro{\ss}e Anzahl von Wiederholungen oder Elementen, um eine gute
Statistik zu bekommen.
\item Benutze die Hilfsfunktion von \code{matlab} (\code{help
commando} oder \code{doc commando}) und das Internet, um
herauszufinden, wie bestimmte \code{matlab} Funktionen zu verwenden
sind und was f\"ur M\"oglichkeiten sie bieten.
Auch zu inhaltlichen Konzepten bietet das Internet oft viele
Antworten!
\end{itemize}
% \begin{itemize}
% \item \"Uberzeuge dich von jeder einzelnen Zeile deines Codes, dass
% sie auch wirklich das macht, was sie machen soll! Teste dies mit
% kleinen Beispielen direkt in der Kommandozeile.
% \item Versuche die L\"osungen der Aufgaben m\"oglichst in
% sinnvolle kleine Funktionen herunterzubrechen.
% Sobald etwas \"ahnliches mehr als einmal berechnet werden soll,
% lohnt es sich eine Funktion daraus zu schreiben!
% \item Teste rechenintensive \code{for} Schleifen, Vektoren, Matrizen
% zuerst mit einer kleinen Anzahl von Wiederholungen oder kleiner
% Gr\"o{\ss}e, und benutze erst am Ende, wenn alles \"uberpr\"uft
% ist, eine gro{\ss}e Anzahl von Wiederholungen oder Elementen, um eine gute
% Statistik zu bekommen.
% \item Benutze die Hilfsfunktion von \code{matlab} (\code{help
% commando} oder \code{doc commando}) und das Internet, um
% herauszufinden, wie bestimmte \code{matlab} Funktionen zu verwenden
% sind und was f\"ur M\"oglichkeiten sie bieten.
% Auch zu inhaltlichen Konzepten bietet das Internet oft viele
% Antworten!
% \end{itemize}
\fi
% \fi

26
likelihood/exercises/mlepropest.m Normal file
View File

@ -0,0 +1,26 @@
m = 2.0; % slope
sigmas = [0.1, 1.0]; % standard deviations
ns = [100, 1000]; % number of data pairs
trials = 1000; % number of data sets
for i = 1:length(sigmas)
sigma = sigmas(i);
for j = 1:length(ns)
n = ns(j);
slopes = zeros(trials, 1);
for k=1:trials
% data pairs:
x = 5.0*rand(n, 1);
y = m*x + sigma*randn(n, 1);
% fit:
slopes(k) = mleslope(x, y);
end
subplot(2, 2, 2*(i-1)+j);
bins = [1.9:0.005:2.1];
hist(slopes, bins);
xlabel('estimated slope');
title(sprintf('sigma=%g, n=%d', sigma, n));
end
end
savefigpdf(gcf, 'mlepropest.pdf', 15, 10);

BIN
likelihood/exercises/mlepropest.pdf Normal file
View File

Binary file not shown.

57
likelihood/exercises/mlestd.m
View File

@ -1,30 +1,35 @@
% draw random numbers:
n = 50;
mu = 3.0;
sigma =2.0;
x = randn(n,1)*sigma+mu;
fprintf(' mean of the data is %.2f\n', mean(x))
fprintf('standard deviation of the data is %.2f\n', std(x))
ns = [50, 1000];
for k = 1:length(ns)
n = ns(k);
% draw random numbers:
x = randn(n,1)*sigma+mu;
fprintf(' mean of the data is %.2f\n', mean(x))
fprintf('standard deviation of the data is %.2f\n', std(x))
% standard deviation as parameter:
psigs = 1.0:0.01:3.0;
% matrix with the probabilities for each x and psigs:
lms = zeros(length(x), length(psigs));
for i=1:length(psigs)
psig = psigs(i);
p = exp(-0.5*((x-mu)/psig).^2.0)/sqrt(2.0*pi)/psig;
lms(:,i) = p;
end
lm = prod(lms, 1); % likelihood
loglm = sum(log(lms), 1); % log likelihood
% standard deviation as parameter:
psigs = 1.0:0.01:3.0;
% matrix with the probabilities for each x and psigs:
lms = zeros(length(x), length(psigs));
for i=1:length(psigs)
psig = psigs(i);
p = exp(-0.5*((x-mu)/psig).^2.0)/sqrt(2.0*pi)/psig;
lms(:,i) = p;
end
lm = prod(lms, 1); % likelihood
loglm = sum(log(lms), 1); % log likelihood
% plot likelihood of standard deviation:
subplot(1, 2, 1);
plot(psigs, lm );
xlabel('standard deviation')
ylabel('likelihood')
subplot(1, 2, 2);
plot(psigs, loglm);
xlabel('standard deviation')
ylabel('log likelihood')
savefigpdf(gcf, 'mlestd.pdf', 15, 5);
% plot likelihood of standard deviation:
subplot(2, 2, 2*k-1);
plot(psigs, lm );
title(sprintf('likelihood n=%d', n));
xlabel('standard deviation')
ylabel('likelihood')
subplot(2, 2, 2*k);
plot(psigs, loglm);
title(sprintf('log-likelihood n=%d', n));
xlabel('standard deviation')
ylabel('log likelihood')
end
savefigpdf(gcf, 'mlestd.pdf', 15, 10);

BIN
likelihood/exercises/mlestd.pdf
View File

Binary file not shown.

52
linearalgebra/code/coordinaterafo2.m
View File

@ -1,52 +0,0 @@
% some vectors:
x = [ -1:0.02:1 ];
y = x*0.5 + 0.1*randn( size(x) );
plot( x, y, '.b' );
hold on
plot( x(10), y(10), '.r' );
% new coordinate system:
%e1 = [ 3/5 4/5 ];
%e2 = [ -4/5 3/5 ];
e1 = [ 3 4 ];
e2 = [ -4 3 ];
me1 = sqrt( e1*e1' );
e1 = e1/me1;
me2 = sqrt( e2*e2' );
e2 = e2/me2;
quiver( 0.0, 0.0 , e1(1), e1(2), 1.0, 'r' )
quiver( 0.0, 0.0 , e2(1), e2(2), 1.0, 'g' )
axis( 'equal' )
% project [x y] onto e1 and e2:
% % the long way:
% nx = zeros( size( x ) ); % new x coordinate
% ny = zeros( size( y ) ); % new y coordinates
% for k=1:length(x)
% xvec = [ x(k) y(k) ];
% nx(k) = xvec * e1';
% ny(k) = xvec * e2';
% end
% plot( nx, ny, '.g' );
% the short way:
%nx = [ x' y' ] * e1';
%nx = [x; y]' * e1';
% nx = e1 * [ x; y ];
% ny = e2 * [ x; y ];
% plot( nx, ny, '.g' );
% even shorter:
n = [e1; e2 ] * [ x; y ];
plot( n(1,:), n(2,:), '.g' );
hold off

4
linearalgebra/code/coordinatetrafo.m
View File

@ -2,8 +2,9 @@
x=[-1:0.02:1]'; % column vector!
y=4*x/3 + 0.1*randn( size( x ) );
plot( x, y, '.b' );
axis( 'equal' );
hold on;
plot( x(1), y(1), '.b' );
axis( 'equal' );
% new coordinate system:
% e1 = [ 3 4 ]
@ -35,3 +36,4 @@ xlabel( 'x' );
ylabel( 'y' );
hold off;
pause;

6
linearalgebra/code/covareigen.m
View File

@ -38,11 +38,11 @@ function covareigen( x, y, pm, w, h )
colormap( 'gray' );
else
% scatter plot:
scatter( x, y, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( x, y ); %, 'b', 'filled', 'MarkerEdgeColor', 'white' );
end
% plot eigenvectors:
quiver( ones( 1, 2 ).*mean( x ), ones( 1, 2 )*mean(y), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
quiver( ones( 1, 2 ).*mean( x ), ones( 1, 2 )*mean(y), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r', 'LineWidth', 3 ); %, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
xlabel( 'x' );
ylabel( 'y' );
@ -76,7 +76,7 @@ function covareigen( x, y, pm, w, h )
contourf( xp, yp, gauss )
end
else
scatter( nc(:,1), nc(:,2), 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( nc(:,1), nc(:,2) ); %, 'b', 'filled', 'MarkerEdgeColor', 'white' );
end
xlabel( 'x' );
ylabel( 'y' );

12
linearalgebra/code/covareigen3.m
View File

@ -15,29 +15,29 @@ function covareigen3( x, y, z )
% scatter plot:
view( 3 );
scatter3( x, y, z, 0.1, 'b', 'filled', 'MarkerEdgeColor', 'blue' );
scatter3( x, y, z, 0.1, 'b', 'filled' ); %, 'MarkerEdgeColor', 'blue' );
xlabel( 'x' );
ylabel( 'y' );
zlabel( 'z' );
grid on;
% plot eigenvectors:
quiver3( ones( 1, 3 ).*mean( x ), ones( 1, 3 )*mean(y), ones( 1, 3 )*mean(z), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', v(3,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
quiver3( ones( 1, 3 ).*mean( x ), ones( 1, 3 )*mean(y), ones( 1, 3 )*mean(z), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', v(3,:).*sqrt(diag(d))', 'r' ); %, 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
%axis( 'equal' );
hold off;
% 2D scatter plots:
subplot( 2, 6, 4 );
scatter( x, y, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( x, y, 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'x' )
ylabel( 'y' )
subplot( 2, 6, 5 );
scatter( x, z, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( x, z, 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'x' )
ylabel( 'z' )
subplot( 2, 6, 6 );
scatter( y, z, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( y, z, 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'y' )
ylabel( 'z' )
@ -51,7 +51,7 @@ function covareigen3( x, y, z )
y = y - mean( y );
% project onto eigenvectors:
nx = [ x y z ] * v(:,inx);
scatter( nx(:,1), nx(:,2), 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( nx(:,1), nx(:,2), 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'ex' )
ylabel( 'ey' )
axis( 'equal' );

4
linearalgebra/code/covareigen3examples.m
View File

@ -1,6 +1,3 @@
scrsz = get( 0, 'ScreenSize' );
set( 0, 'DefaultFigurePosition', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ] );
n = 10000;
% three distributions:
@ -17,3 +14,4 @@ for k = 1:4
end
f = figure( 1 );
covareigen3( x, y, z );
pause;

9
linearalgebra/code/covareigenexamples.m
View File

@ -1,6 +1,3 @@
scrsz = get( 0, 'ScreenSize' );
set( 0, 'DefaultFigurePosition', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ] );
% correlation coefficients:
n = 10000;
x = randn( n, 1 );
@ -10,9 +7,9 @@ for r = 0.01:0.19:1
clf( f );
y = r*x + sqrt(1-r^2)*randn( n, 1 );
covareigen( x, y, 0, 5.0, 3.0 );
key = waitforbuttonpress;
%key = waitforbuttonpress;
pause( 1.0 );
end
return
% two distributions:
n = 10000;
@ -21,14 +18,12 @@ y1 = randn( n/2, 1 );
x2 = randn( n/2, 1 );
y2 = randn( n/2, 1 );
f = figure( 1 );
pause( 'on' );
for d = 0:1:5
fprintf( 'Distance = %g\n', d );
clf( f );
d2 = d / sqrt( 2.0 );
x = [ x1; x2 ];
y = [ y1+d2; y2-d2 ];
scrsz = get(0,'ScreenSize');
covareigen( x, y, 0, 10.0, 7.0 );
%key = waitforbuttonpress;
pause( 1.0 );

28
linearalgebra/code/covarmatrix.m Normal file
View File

@ -0,0 +1,28 @@
% covariance in 2D:
x = 2.0*randn(1000, 1);
y = 0.5*x + 0.2*randn(1000, 1);
scatter(x, y)
var(x)
var(y)
cov( [x, y] )
pause(1.0)
% covariance of independent coordinates in 2D:
x = 2.0*randn(1000, 1);
y = 0.5*randn(1000, 1);
scatter(x, y)
var(x)
var(y)
cov( [x, y] )
pause(1.0)
% covariance in 3D:
x = 2.0*randn(1000, 1);
y = -0.5*x + 0.2*randn(1000, 1);
z = 2.0*x + 0.6*y + 2.0*randn(1000, 1);
scatter3(x, y, z, 'filled')
var(x)
var(y)
var(z)
cov( [x, y, z] )
pause()

36
linearalgebra/code/covarmatrixeigen.m Normal file
View File

@ -0,0 +1,36 @@
% covariance in 2D:
x = 2.0*randn(1000, 1);
y = 0.5*x + 0.6*randn(1000, 1);
scatter(x, y)
cv = cov( [x, y] );
[v d] = eig(cv);
hold on
quiver([0 0], [0 0], v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r', 'LineWidth', 2);
axis( 'equal' );
hold off
pause
% covariance of independent coordinates in 2D:
x = 1.5*randn(1000, 1);
y = 0.8*randn(1000, 1);
scatter(x, y)
cv = cov( [x, y] );
[v d] = eig(cv);
hold on
quiver([0 0], [0 0], v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r', 'LineWidth', 2);
axis( 'equal' );
hold off
pause
% covariance in 3D:
x = 2.0*randn(1000, 1);
y = -0.5*x + 0.2*randn(1000, 1);
z = 2.0*x + 0.6*y + 2.0*randn(1000, 1);
scatter3(x, y, z)
cv = cov( [x, y, z] );
[v d] = eig(cv);
hold on
quiver3([0 0 0], [0 0 0], [0 0 0], v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', v(3,:).*sqrt(diag(d))', 'r', 'LineWidth', 2);
axis( 'equal' );
hold off
pause()

2
linearalgebra/code/matrixbox.m
View File

@ -41,5 +41,5 @@ function matrixbox( m, s )
text( 0.8, 0.1, sprintf( '%.3g', m(2,1) ), 'Units', 'normalized' )
text( 0.9, 0.1, sprintf( '%.3g', m(2,2) ), 'Units', 'normalized' )
title( s );
waitforbuttonpress;
pause(1.0);
end

6
linearalgebra/code/pca2d.m
View File

@ -39,11 +39,11 @@ function pca2d( x, y, pm, w, h )
colormap( 'gray' );
else
% scatter plot:
scatter( x, y, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( x, y, 'b', 'filled' );%, 'MarkerEdgeColor', 'white' );
end
% plot eigenvectors:
quiver( ones( 1, 2 ).*mean( x ), ones( 1, 2 )*mean(y), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
quiver( ones( 1, 2 ).*mean( x ), ones( 1, 2 )*mean(y), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r' ); %, 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
xlabel( 'x' );
ylabel( 'y' );
@ -80,7 +80,7 @@ function pca2d( x, y, pm, w, h )
end
else
% scatter plot:
scatter( nc(:,1), nc(:,2), 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( nc(:,1), nc(:,2), 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
end
xlabel( 'n_x' );
ylabel( 'n_y' );

14
linearalgebra/code/pca2dexamples.m
View File

@ -1,6 +1,3 @@
scrsz = get( 0, 'ScreenSize' );
set( 0, 'DefaultFigurePosition', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ] );
% correlation coefficients:
n = 10000;
x = randn( n, 1 );
@ -10,7 +7,7 @@ for r = 0.01:0.19:1
clf( f );
y = r*x + sqrt(1-r^2)*randn( n, 1 );
pca2d( x, y, 0, 5.0, 3.0 );
waitforbuttonpress;
pause( 1.0 );
end
% two distributions:
@ -20,14 +17,11 @@ y1 = randn( n/2, 1 );
x2 = randn( n/2, 1 );
y2 = randn( n/2, 1 );
f = figure( 1 );
pause( 'on' );
for d = 0:1:5
fprintf( 'Distance = %g\n', d );
clf( f );
d2 = d / sqrt( 2.0 );
x = [ x1; x2 ];
y = [ y1+d2; y2-d2 ];
scrsz = get(0,'ScreenSize');
x = [ x1+d; x2-d ];
y = [ y1+d; y2-d ];
pca2d( x, y, 0, 10.0, 7.0 );
waitforbuttonpress;
pause( 1.0 );
end

13
linearalgebra/code/pca3d.m
View File

@ -15,29 +15,29 @@ function pca3d( x, y, z )
% scatter plot:
view( 3 );
scatter3( x, y, z, 0.1, 'b', 'filled', 'MarkerEdgeColor', 'blue' );
scatter3( x, y, z ); % , 1.0, 'b', 'filled' ); %, 'MarkerEdgeColor', 'blue' );
xlabel( 'x' );
ylabel( 'y' );
zlabel( 'z' );
grid on;
% plot eigenvectors:
quiver3( ones( 1, 3 ).*mean( x ), ones( 1, 3 )*mean(y), ones( 1, 3 )*mean(z), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', v(3,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
quiver3( ones( 1, 3 ).*mean( x ), ones( 1, 3 )*mean(y), ones( 1, 3 )*mean(z), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', v(3,:).*sqrt(diag(d))', 'r' ); %, 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
%axis( 'equal' );
hold off;
% 2D scatter plots:
subplot( 2, 6, 4 );
scatter( x, y, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( x, y, 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'x' )
ylabel( 'y' )
subplot( 2, 6, 5 );
scatter( x, z, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( x, z, 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'x' )
ylabel( 'z' )
subplot( 2, 6, 6 );
scatter( y, z, 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( y, z, 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'y' )
ylabel( 'z' )
@ -51,10 +51,11 @@ function pca3d( x, y, z )
y = y - mean( y );
% project onto eigenvectors:
nx = [ x y z ] * v(:,inx);
scatter( nx(:,1), nx(:,2), 'b', 'filled', 'MarkerEdgeColor', 'white' );
scatter( nx(:,1), nx(:,2), 'b', 'filled' ); %, 'MarkerEdgeColor', 'white' );
xlabel( 'ex' )
ylabel( 'ey' )
axis( 'equal' );
hold off;
end

4
linearalgebra/code/pca3dexamples.m
View File

@ -1,6 +1,3 @@
scrsz = get( 0, 'ScreenSize' );
set( 0, 'DefaultFigurePosition', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ] );
n = 10000;
% three distributions:
@ -17,3 +14,4 @@ for k = 1:4
end
f = figure( 1 );
pca3d( x, y, z );
pause

2
linearalgebra/code/simplematrixbox.m
View File

@ -17,5 +17,5 @@ function simplematrixbox( a, s )
xlim( [-2 2 ] )
ylim( [-2 2] )
title( s );
waitforbuttonpress;
pause(1.0);
end

40
linearalgebra/code/spikesortingwave.m
View File

@ -5,21 +5,22 @@ load( 'extdata' );
dt = time( 2) - time(1);
tinx = round(spiketimes/dt)+1;
% plot voltage trace with dettected spikes:
% plot voltage trace with detected spikes:
figure( 1 );
clf;
plot( time, voltage, '-b' )
hold on
scatter( time(tinx), voltage(tinx), 'r', 'filled' );
xlabel( 'time [ms]' );
xlabel( 'time [s]' );
ylabel( 'voltage' );
xlim([0.1, 0.4])
hold off
% spike waveform snippets:
w = ceil( 0.005/dt );
vs = [];
vs = zeros(length(tinx), 2*w);
for k=1:length(tinx)
vs = [ vs; voltage(tinx(k)-w:tinx(k)+w-1) ];
vs(k,:) = voltage(tinx(k)-w:tinx(k)+w-1);
end
ts = time(1:size(vs,2));
ts = ts - ts(floor(length(ts)/2));
@ -78,9 +79,9 @@ kx = ones( size( nx ) );
kx(nx<kthresh) = 2;
subplot( 4, 1, 2 );
scatter( nx(kx==1), ny(kx==1), 'r', 'filled', 'MarkerEdgeColor', 'white' );
scatter( nx(kx==1), ny(kx==1), 10.0, 'r', 'filled' ); %, 'MarkerEdgeColor', 'white' );
hold on;
scatter( nx(kx==2), ny(kx==2), 'g', 'filled', 'MarkerEdgeColor', 'white' );
scatter( nx(kx==2), ny(kx==2), 10.0, 'g', 'filled' ); %, 'MarkerEdgeColor', 'white' );
hold off;
xlabel( 'projection onto eigenvector 1' );
ylabel( 'projection onto eigenvector 2' );
@ -113,6 +114,29 @@ hold off
% spike trains:
figure( 1 );
hold on;
scatter( time(tinx(kinx1)), voltage(tinx(kinx1)), 100.0, 'r', 'filled' );
scatter( time(tinx(kinx2)), voltage(tinx(kinx2)), 100.0, 'g', 'filled' );
scatter( time(tinx(kinx1)), voltage(tinx(kinx1)), 10.0, 'r', 'filled' );
scatter( time(tinx(kinx2)), voltage(tinx(kinx2)), 10.0, 'g', 'filled' );
hold off;
% ISIs:
figure(4);
clf;
subplot(1, 3, 1)
allisis = diff(spiketimes);
hist(allisis, 20);
title('all spikes')
xlabel('ISI [ms]')
subplot(1, 3, 2)
spikes1 = time(tinx(kinx1));
isis1 = diff(spikes1);
hist(isis1, 20);
title('neuron 1')
xlabel('ISI [ms]')
subplot(1, 3, 3)
spikes2 = time(tinx(kinx2));
isis2 = diff(spikes2);
hist(isis2, 20);
title('neuron 2')
xlabel('ISI [ms]')
pause

34
linearalgebra/exercises/Makefile Normal file
View File

@ -0,0 +1,34 @@
TEXFILES=$(wildcard exercises??.tex)
EXERCISES=$(TEXFILES:.tex=.pdf)
SOLUTIONS=$(EXERCISES:exercises%=solutions%)
.PHONY: pdf exercises solutions watch watchexercises watchsolutions clean
pdf : $(SOLUTIONS) $(EXERCISES)
exercises : $(EXERCISES)
solutions : $(SOLUTIONS)
$(SOLUTIONS) : solutions%.pdf : exercises%.tex instructions.tex
{ echo "\\documentclass[answers,12pt,a4paper,pdftex]{exam}"; sed -e '1d' $<; } > $(patsubst %.pdf,%.tex,$@)
pdflatex -interaction=scrollmode $(patsubst %.pdf,%.tex,$@) | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $(patsubst %.pdf,%.tex,$@) || true
rm $(patsubst %.pdf,%,$@).[!p]*
$(EXERCISES) : %.pdf : %.tex instructions.tex
pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true
watch :
while true; do ! make -q pdf && make pdf; sleep 0.5; done
watchexercises :
while true; do ! make -q exercises && make exercises; sleep 0.5; done
watchsolutions :
while true; do ! make -q solutions && make solutions; sleep 0.5; done
clean :
rm -f *~ *.aux *.log *.out
cleanup : clean
rm -f $(SOLUTIONS) $(EXERCISES)

BIN
linearalgebra/exercises/UT_WBMW_Black_RGB.pdf Normal file
View File

Binary file not shown.

19
linearalgebra/exercises/covariance.m Normal file
View File

@ -0,0 +1,19 @@
n = 1000;
x = 2.0*randn(n, 1);
z = 3.0*randn(n, 1);
rs = [-1.0:0.25:1.0];
for k = 1:length(rs)
r = rs(k);
y = r*x + sqrt(1.0-r^2)*z;
r
cv = cov([x y])
corrcoef([x y])
subplot(3, 3, k)
scatter(x, y, 'filled', 'MarkerEdgeColor', 'white')
title(sprintf('r=%g', r))
xlim([-10, 10])
ylim([-20, 20])
text(-8, -15, sprintf('cov(x,y)=%.2f', cv(1, 2)))
end
savefigpdf(gcf, 'covariance.pdf', 20, 20);

BIN
linearalgebra/exercises/covariance.pdf Normal file
View File

Binary file not shown.

244
linearalgebra/exercises/exercises01.tex Normal file
View File

@ -0,0 +1,244 @@
\documentclass[12pt,a4paper,pdftex]{exam}
\usepackage[english]{babel}
\usepackage{pslatex}
\usepackage[mediumspace,mediumqspace,Gray]{SIunits} % \ohm, \micro
\usepackage{xcolor}
\usepackage{graphicx}
\usepackage[breaklinks=true,bookmarks=true,bookmarksopen=true,pdfpagemode=UseNone,pdfstartview=FitH,colorlinks=true,citecolor=blue]{hyperref}
%%%%% layout %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[left=20mm,right=20mm,top=25mm,bottom=25mm]{geometry}
\pagestyle{headandfoot}
\ifprintanswers
\newcommand{\stitle}{: Solutions}
\else
\newcommand{\stitle}{}
\fi
\header{{\bfseries\large Exercise\stitle}}{{\bfseries\large PCA}}{{\bfseries\large January 7th, 2019}}
\firstpagefooter{Prof. Dr. Jan Benda}{Phone: 29 74573}{Email:
jan.benda@uni-tuebingen.de}
\runningfooter{}{\thepage}{}
\setlength{\baselineskip}{15pt}
\setlength{\parindent}{0.0cm}
\setlength{\parskip}{0.3cm}
\renewcommand{\baselinestretch}{1.15}
%%%%% listings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{listings}
\lstset{
language=Matlab,
basicstyle=\ttfamily\footnotesize,
numbers=left,
numberstyle=\tiny,
title=\lstname,
showstringspaces=false,
commentstyle=\itshape\color{darkgray},
breaklines=true,
breakautoindent=true,
columns=flexible,
frame=single,
xleftmargin=1em,
xrightmargin=1em,
aboveskip=10pt
}
%%%%% math stuff: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{bm}
\usepackage{dsfont}
\newcommand{\naZ}{\mathds{N}}
\newcommand{\gaZ}{\mathds{Z}}
\newcommand{\raZ}{\mathds{Q}}
\newcommand{\reZ}{\mathds{R}}
\newcommand{\reZp}{\mathds{R^+}}
\newcommand{\reZpN}{\mathds{R^+_0}}
\newcommand{\koZ}{\mathds{C}}
%%%%% page breaks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\continue}{\ifprintanswers%
\else
\vfill\hspace*{\fill}$\rightarrow$\newpage%
\fi}
\newcommand{\continuepage}{\ifprintanswers%
\newpage
\else
\vfill\hspace*{\fill}$\rightarrow$\newpage%
\fi}
\newcommand{\newsolutionpage}{\ifprintanswers%
\newpage%
\else
\fi}
%%%%% new commands %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\qt}[1]{\textbf{#1}\\}
\newcommand{\pref}[1]{(\ref{#1})}
\newcommand{\extra}{--- Zusatzaufgabe ---\ \mbox{}}
\newcommand{\code}[1]{\texttt{#1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\input{instructions}
\begin{questions}
\question \qt{Covariance and correlation coefficient\vspace{-3ex}}
\begin{parts}
\part Generate two vectors $x$ and $z$ with $n=1000$ Gausian distributed random numbers.
\part Compute $y$ as a linear combination of $x$ and $z$ according to
\[ y = r \cdot x + \sqrt{1-r^2}\cdot z \]
where $r$ is a parameter $-1 \le r \le 1$.
What does $r$ do?
\part Plot a scatter plot of $y$ versus $x$ for about 10 different values of $r$.
What do you observe?
\part Also compute the covariance matrix and the correlation
coefficient matrix between $x$ and $y$ (functions \texttt{cov} and
\texttt{corrcoef}). How do these matrices look like for different
values of $r$? How do the values of the matrices change if you generate
$x$ and $z$ with larger variances?
\begin{solution}
\lstinputlisting{covariance.m}
\includegraphics[width=0.8\textwidth]{covariance}
\end{solution}
\part Do the same analysis (scatter plot, covariance, and correlation coefficient)
for \[ y = x^2 + 0.5 \cdot z \]
Are $x$ and $y$ really independent?
\begin{solution}
\lstinputlisting{nonlinearcorrelation.m}
\includegraphics[width=0.8\textwidth]{nonlinearcorrelation}
\end{solution}
\end{parts}
\question \qt{Principal component analysis in 2D\vspace{-3ex}}
\begin{parts}
\part Generate $n=1000$ pairs $(x,y)$ of Gaussian distributed random numbers such
that all $x$ values have zero mean, half of the $y$ values have mean $+d$
and the other half mean $-d$, with $d \ge0$.
\part Plot scatter plots of the pairs $(x,y)$ for $d=0$, 1, 2, 3, 4 and 5.
Also plot a histogram of the $x$ values.
\part Apply PCA on the data and plot a histogram of the data projected onto
the PCA axis with the largest eigenvalue.
What do you observe?
\begin{solution}
\lstinputlisting{pca2d.m}
\includegraphics[width=0.8\textwidth]{pca2d-2}
\end{solution}
\end{parts}
\newsolutionpage
\question \qt{Principal component analysis in 3D\vspace{-3ex}}
\begin{parts}
\part Generate $n=1000$ triplets $(x,y,z)$ of Gaussian distributed random numbers such
that all $x$ values have zero mean, half of the $y$ and $z$ values have mean $+d$
and the other half mean $-d$, with $d \ge0$.
\part Plot 3D scatter plots of the pairs $(x,y)$ for $d=0$, 1, 2, 3, 4 and 5.
Also plot a histogram of the $x$ values.
\part Apply PCA on the data and plot a histogram of the data projected onto
the PCA axis with the largest eigenvalue.
What do you observe?
\begin{solution}
\lstinputlisting{pca3d.m}
\includegraphics[width=0.8\textwidth]{pca3d-2}
\end{solution}
\end{parts}
\continuepage
\question \qt{Spike sorting}
Extracellular recordings often pick up action potentials originating
from more than a single neuron. In case the waveforms of the action
potentials differ between the neurons one could assign each action
potential to the neuron it originated from. This process is called
``spike sorting''. Here we explore this method on a simulated
recording that contains action potentials from two different
neurons.
\begin{parts}
\part Load the data from the file \texttt{extdata.mat}. This file
contains the voltage trace of the recording (\texttt{voltage}),
the corresponding time vector in seconds (\texttt{time}), and a
vector containing the times of the peaks of detected action
potentials (\texttt{spiketimes}). Further, and in contrast to real
data, the waveforms of the actionpotentials of the two neurons
(\texttt{waveform1} and \texttt{waveform2}) and the corresponding
time vector (\texttt{waveformt}) are also contained in the file.
\part Plot the voltage trace and mark the peaks of the detected
action potentials (using \texttt{spiketimes}). Zoom into the plot
and look whether you can differentiate between two different
waveforms of action potentials. How do they differ?
\part Cut out the waveform of each action potential (5\,ms before
and after the peak). Plot all these snippets in a single
plot. Can you differentiate the two actionpotential waveforms?
\begin{solution}
\mbox{}\\[-3ex]\hspace*{5em}
\includegraphics[width=0.8\textwidth]{spikesorting1}
\end{solution}
\newsolutionpage
\part Apply PCA on the waveform snippets. That is compute the
eigenvalues and eigenvectors of their covariance matrix, which is
a $n \times n$ matrix, with $n$ being the number of data points
contained in a single waveform snippet. Plot the sorted
eigenvalues (the ``eigenvalue spectrum''). How many eigenvalues
are clearly larger than zero?
\begin{solution}
\mbox{}\\[-3ex]\hspace*{5em}
\includegraphics[width=0.8\textwidth]{spikesorting2}
\end{solution}
\part Plot the two eigenvectors (``features'') with the two
largest eigenvalues as a function of time.
\begin{solution}
\mbox{}\\[-3ex]\hspace*{5em}
\includegraphics[width=0.8\textwidth]{spikesorting3}
\end{solution}
\part Project the waveform snippets onto these two eigenvectors
and display them with a scatter plot. What do you observe? Can you
imagine how to separate two ``clouds'' of data points
(``clusters'')?
\newsolutionpage
\part Think about a very simply way how to separate the two
clusters. Generate a vector whose elements label the action
potentials, e.g. that contains '1' for all snippets belonging to
the one cluster and '2' for the waveforms of the other
cluster. Use this vector to mark the two clusters in the previous
plot with two different colors.
\begin{solution}
\mbox{}\\[-3ex]\hspace*{5em}
\includegraphics[width=0.8\textwidth]{spikesorting4}
\end{solution}
\part Plot the waveform snippets of each cluster together with the
true waveform obtained from the data file. Do they match?
\begin{solution}
\mbox{}\\[-3ex]\hspace*{5em}
\includegraphics[width=0.8\textwidth]{spikesorting5}
\end{solution}
\newsolutionpage
\part Mark the action potentials in the recording according to
their cluster identity.
\begin{solution}
\mbox{}\\[-3ex]\hspace*{5em}
\includegraphics[width=0.8\textwidth]{spikesorting6}
\end{solution}
\part Compute interspike-interval histograms of all the (unsorted)
action potentials, and of each of the two neurons. What do they
tell you?
\begin{solution}
\mbox{}\\[-3ex]\hspace*{5em}
\includegraphics[width=0.8\textwidth]{spikesorting7}
\lstinputlisting{spikesorting.m}
\end{solution}
\end{parts}
\end{questions}
\end{document}

BIN
linearalgebra/exercises/extdata.mat Normal file
View File

Binary file not shown.

6
linearalgebra/exercises/instructions.tex Normal file
View File

@ -0,0 +1,6 @@
\vspace*{-7.8ex}
\begin{center}
\textbf{\Large Introduction to Scientific Computing}\\[2.3ex]
{\large Jan Grewe, Jan Benda}\\[-3ex]
Neuroethology Lab \hfill --- \hfill Institute for Neurobiology \hfill --- \hfill \includegraphics[width=0.28\textwidth]{UT_WBMW_Black_RGB} \\
\end{center}

9
linearalgebra/exercises/nonlinearcorrelation.m Normal file
View File

@ -0,0 +1,9 @@
n = 1000;
x = randn(n, 1);
z = randn(n, 1);
y = x.^2 + 0.5*z;
scatter(x, y)
cov([x y])
r = corrcoef([x y])
text(-2, 8, sprintf('r=%.2f', r(1, 2)))
savefigpdf(gcf, 'nonlinearcorrelation.pdf', 15, 10);

BIN
linearalgebra/exercises/nonlinearcorrelation.pdf Normal file
View File

Binary file not shown.

BIN
linearalgebra/exercises/pca2d-2.pdf Normal file
View File

Binary file not shown.

39
linearalgebra/exercises/pca2d.m Normal file
View File

@ -0,0 +1,39 @@
n = 1000;
ds = [0:5];
for k = 1:length(ds)
d = ds(k);
% generate data:
x = randn(n, 1);
y = randn(n, 1);
y(1:n/2) = y(1:n/2) - d;
y(1+n/2:end) = y(1+n/2:end) + d;
% scatter plot of data:
subplot(2, 2, 1)
scatter(x, y, 'filled', 'MarkerEdgeColor', 'white')
title(sprintf('d=%g', d))
xlabel('x')
ylabel('y')
% histogram of data:
subplot(2, 2, 3)
hist(x, 20)
xlabel('x')
% pca:
cv = cov([x y]);
[ev en] = eig(cv);
[en inx] = sort(diag(en), 'descend');
nc = [x y]*ev(:,inx);
% scatter in new coordinates:
subplot(2, 2, 2)
scatter(nc(:, 1), nc(:, 2), 'filled', 'MarkerEdgeColor', 'white')
title(sprintf('d=%g', d))
xlabel('pca 1')
ylabel('pca 2')
% histogram of data:
subplot(2, 2, 4)
hist(nc(:, 1), 20)
xlabel('pca 1')
if d == 2
savefigpdf(gcf, 'pca2d-2.pdf', 15, 15);
end
pause(1.0)
end

BIN
linearalgebra/exercises/pca3d-2.pdf Normal file
View File

Binary file not shown.

44
linearalgebra/exercises/pca3d.m Normal file
View File

@ -0,0 +1,44 @@
n = 1000;
ds = [0:5];
for k = 1:length(ds)
d = ds(k);
% generate data:
x = randn(n, 1);
y = randn(n, 1);
z = randn(n, 1);
y(1:n/2) = y(1:n/2) - d;
y(1+n/2:end) = y(1+n/2:end) + d;
z(1:n/2) = z(1:n/2) - d;
z(1+n/2:end) = z(1+n/2:end) + d;
% scatter plot of data:
subplot(2, 2, 1)
scatter3(x, y, z, 'filled', 'MarkerEdgeColor', 'white')
title(sprintf('d=%g', d))
xlabel('x')
ylabel('y')
zlabel('z')
% histogram of data:
subplot(2, 2, 3)
hist(x, 20)
xlabel('x')
% pca:
cv = cov([x y z]);
[ev en] = eig(cv);
[en inx] = sort(diag(en), 'descend');
nc = [x y z]*ev(:,inx);
% scatter in new coordinates:
subplot(2, 2, 2)
scatter3(nc(:, 1), nc(:, 2), nc(:, 3), 'filled', 'MarkerEdgeColor', 'white')
title(sprintf('d=%g', d))
xlabel('pca 1')
ylabel('pca 2')
zlabel('pca 3')
% histogram of data:
subplot(2, 2, 4)
hist(nc(:, 1), 20)
xlabel('pca 1')
if d == 2
savefigpdf(gcf, 'pca3d-2.pdf', 15, 15);
end
pause(1.0)
end

28
linearalgebra/exercises/savefigpdf.m Normal file
View File

@ -0,0 +1,28 @@
function savefigpdf(fig, name, width, height)
% Saves figure fig in pdf file name.pdf with appropriately set page size
% and fonts
% default width:
if nargin < 3
width = 11.7;
end
% default height:
if nargin < 4
height = 9.0;
end
% paper:
set(fig, 'PaperUnits', 'centimeters');
set(fig, 'PaperSize', [width height]);
set(fig, 'PaperPosition', [0.0 0.0 width height]);
set(fig, 'Color', 'white')
% font:
set(findall(fig, 'type', 'axes'), 'FontSize', 12)
set(findall(fig, 'type', 'text'), 'FontSize', 12)
% save:
saveas(fig, name, 'pdf')
end

161
linearalgebra/exercises/spikesorting.m Normal file
View File

@ -0,0 +1,161 @@
% load data into time, voltage and spiketimes
x = load('extdata');
time = x.time;
voltage = x.voltage;
spiketimes = x.spiketimes;
waveformt = x.waveformt;
waveform1 = x.waveform1;
waveform2 = x.waveform2;
% indices into voltage trace of spike times:
dt = time(2) - time(1);
tinx = round(spiketimes/dt) + 1;
% plot voltage trace with detected spikes:
figure(1);
plot(time, voltage, '-b')
hold on
scatter(time(tinx), voltage(tinx), 'r', 'filled');
xlabel('time [s]');
ylabel('voltage');
xlim([0.1, 0.4]) % zoom in
hold off
% spike waveform snippets:
snippetwindow = 0.005; % milliseconds
w = ceil(snippetwindow/dt);
vs = zeros(length(tinx), 2*w);
for k=1:length(tinx)
vs(k,:) = voltage(tinx(k)-w:tinx(k)+w-1);
end
% time axis for snippets:
ts = time(1:size(vs,2));
ts = ts - ts(floor(length(ts)/2));
% plot all snippets:
figure(2);
hold on
plot(1000.0*ts, vs, '-b');
title('spike snippets')
xlabel('time [ms]');
ylabel('voltage');
hold off
savefigpdf(gcf, 'spikesorting1.pdf', 12, 6);
% pca:
cv = cov(vs);
[ev , ed] = eig(cv);
[d, dinx] = sort(diag(ed), 'descend');
% features:
figure(2)
subplot(1, 2, 1);
plot(1000.0*ts, ev(:,dinx(1)), 'r', 'LineWidth', 2);
xlabel('time [ms]');
ylabel('eigenvector 1');
subplot(1, 2, 2);
plot(1000.0*ts, ev(:,dinx(2)), 'g', 'LineWidth', 2);
xlabel('time [ms]');
ylabel('eigenvector 2');
savefigpdf(gcf, 'spikesorting2.pdf', 12, 6);
% plot covariance matrix:
figure(3);
subplot(1, 2, 1);
imagesc(cv);
xlabel('time bin');
ylabel('time bin');
title('covariance matrix');
caxis([-0.1 0.1])
% spectrum of eigenvalues:
subplot(1, 2, 2);
scatter(1:length(d), d, 'b', 'filled');
title('spectrum');
xlabel('index');
ylabel('eigenvalue');
savefigpdf(gcf, 'spikesorting3.pdf', 12, 6);
% project onto eigenvectors:
nx = vs * ev(:,dinx(1));
ny = vs * ev(:,dinx(2));
% clustering (two clusters):
%kx = kmeans([nx, ny], 2);
% nx smaller or greater a threshold:
kthresh = 1.6;
kx = ones(size(nx));
kx(nx<kthresh) = 2;
% plot pca coordinates:
figure(4)
scatter(nx(kx==1), ny(kx==1), 'r', 'filled');
hold on;
scatter(nx(kx==2), ny(kx==2), 'g', 'filled');
hold off;
xlabel('projection onto eigenvector 1');
ylabel('projection onto eigenvector 2');
savefigpdf(gcf, 'spikesorting4.pdf', 12, 10);
% show sorted spike waveforms:
figure(5)
subplot(1, 2, 1);
hold on
plot(1000.0*ts, vs(kx==1,:), '-r');
plot(1000.0*waveformt, waveform2, '-k', 'LineWidth', 2);
xlim([1000.0*ts(1) 1000.0*ts(end)])
xlabel('time [ms]');
ylabel('waveform 1');
hold off
subplot(1, 2, 2);
hold on
plot(1000.0*ts, vs(kx==2,:), '-g');
plot(1000.0*waveformt, waveform1, '-k', 'LineWidth', 2);
xlim([1000.0*ts(1) 1000.0*ts(end)])
xlabel('time [ms]');
ylabel('waveform 2');
hold off
savefigpdf(gcf, 'spikesorting5.pdf', 12, 6);
% mark neurons in recording:
figure(1);
hold on;
scatter(time(tinx(kinx1)), voltage(tinx(kinx1)), 'r', 'filled');
scatter(time(tinx(kinx2)), voltage(tinx(kinx2)), 'g', 'filled');
hold off;
savefigpdf(gcf, 'spikesorting6.pdf', 12, 6);
% ISIs:
figure(7);
subplot(1, 3, 1)
allisis = diff(spiketimes);
bins = [0:0.005:0.2];
[h, b] = hist(allisis, bins);
bar(1000.0*b, h/sum(h)/mean(diff(b)))
title('all spikes')
xlabel('ISI [ms]')
xlim([0, 200.0])
subplot(1, 3, 2)
spikes1 = time(tinx(kx==1));
isis1 = diff(spikes1);
[h, b] = hist(isis1, bins);
bar(1000.0*b, h/sum(h)/mean(diff(b)))
title('neuron 1')
xlabel('ISI [ms]')
xlim([0, 200.0])
subplot(1, 3, 3)
spikes2 = time(tinx(kx==2));
isis2 = diff(spikes2);
[h, b] = hist(isis2, bins);
bar(1000.0*b, h/sum(h)/mean(diff(b)))
title('neuron 2')
xlabel('ISI [ms]')
xlim([0, 200.0])
savefigpdf(gcf, 'spikesorting7.pdf', 12, 6);

BIN
linearalgebra/exercises/spikesorting1.pdf Normal file
View File

Binary file not shown.

BIN
linearalgebra/exercises/spikesorting2.pdf Normal file
View File

Binary file not shown.

BIN
linearalgebra/exercises/spikesorting3.pdf Normal file
View File

Binary file not shown.

BIN
linearalgebra/exercises/spikesorting4.pdf Normal file
View File

Binary file not shown.

BIN
linearalgebra/exercises/spikesorting5.pdf Normal file
View File

Binary file not shown.

BIN
linearalgebra/exercises/spikesorting6.pdf Normal file
View File

Binary file not shown.

BIN
linearalgebra/exercises/spikesorting7.pdf Normal file
View File

Binary file not shown.

60
projects/README
View File

@ -1,70 +1,92 @@
project_adaptation_fit
OK
OK, medium
Add plotting of cost function
project_eod
Needs to be checked!
OK, medium - difficult
project_eyetracker
OK
OK, difficult
no statistics, but kmeans
project_fano_slope
OK
OK, difficult
project_fano_test
OK
OK -
project_fano_time
OK
OK, medium-difficult
project_ficurves
OK
OK, medium
Maybe add correlation test or fit statistics
project_input_resistance
What is the problem with this project?
medium
What is the problem with this project? --> No difference between segments
Improve questions
project_isicorrelations
Need to program a solution!
medium-difficult
Need to finish solution
project_isipdffit
Too technical
project_lif
OK
OK, difficult
no statistics
project_mutualinfo
OK
OK, medium
project_noiseficurves
OK
OK, simple-medium
no statistics
project_numbers
simple
We might add some more involved statistical analysis
project_pca_natural_images
Needs PCA...
medium
Make a solution (->Lukas)
project_photoreceptor
OK - text needs to be improved
Maybe also add how responses are influenced by unstable resting potential
Maybe more cells...
OK, simple
project_populationvector
difficult
OK
project_qvalues
-
Interesting! But needs solution.
project_random_walk
simple-medium
Improve it! Provide code exmaples for plotting world and making movies
project_serialcorrelation
OK
OK, simple-medium
project_spectra
-
Needs improvements and a solution
project_stimulus_reconstruction
OK Fix equation?
OK, difficult
Add specific hints for statistics
project_vector_strength
OK Maybe add something for explainiong the vector strength (average unit vector).
OK, medium-difficult
Maybe add something for explaining the vector strength (average unit vector).
Check text (d)
Add statisitcs for (e)
Peter power:
medium
Marius monkey data:
medium-difficult

6
projects/evaluation.tex
View File

@ -17,14 +17,14 @@
\begin{document}
\sffamily
\section*{Scientific computing WS17/18}
\section*{Scientific computing WS18/19}
\noindent
\begin{tabular}{|p{0.16\textwidth}|p{0.1\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|p{0.024\textwidth}|l|l|}
\hline
Name & Project & \multicolumn{5}{l|}{Project} & \multicolumn{3}{l|}{Figures} & \multicolumn{8}{l|}{Code} & Sum & Grade \rule{0.pt}{3ex} \\
& & \rot{Understanding} & \rot{Presentation}\rot{Structure, didactics} & \rot{Completed} & \rot{Bonus:}\rot{Context} & \rot{Bonus:}\rot{Own ideas} & \rot{Readability}\rot{font size, etc.} & \rot{Labels}\rot{and units} & \rot{Design} & \rot{Runs} & \rot{Names of}\rot{vars and funcs} & \rot{Style: white space}\rot{Spaghetti, structure} & \rot{Functions: used,}\rot{purpose, standalone} & \rot{Docu script} & \rot{Docu functions} & \rot{Figures saved} & \rot{Bonus:}\rot{Interesting code} & & \\
& & 0--3 & 0--2 & 0--2 & 0, 1 & 0--2 & 0--2 & 0--2 & 0--2 & 0--2 & 0--2 & 0--2 & 0--3 & 0--2 & 0--2 & 0, 2 & 0--2 & 28 & \\ \hline
& & \rot{Understanding} & \rot{Presentation}\rot{Structure, didactics} & \rot{Completed} & \rot{Bonus:}\rot{Context} & \rot{Bonus:}\rot{Own ideas} & \rot{Readability}\rot{font size, etc.} & \rot{Labels, units}\rot{and statistical measures} & \rot{Design} & \rot{Runs} & \rot{Names of}\rot{vars and funcs} & \rot{Style: white space}\rot{Spaghetti, structure} & \rot{Functions: used,}\rot{purpose, standalone} & \rot{Docu script} & \rot{Docu functions} & \rot{Figures saved} & \rot{Bonus:}\rot{Interesting code} & & \\
& & 0--6 & 0--2 & 0--2 & 0, 1 & 0--2 & 0--2 & 0--2 & 0--2 & 0--2 & 0--2 & 0--2 & 0--3 & 0--2 & 0--2 & 0, 2 & 0--2 & 28 & \\ \hline
\num & & & & & & & & & & & & & & & & & & & \\ \hline
\num & & & & & & & & & & & & & & & & & & & \\ \hline
\num & & & & & & & & & & & & & & & & & & & \\ \hline

28
projects/instructions.tex
View File

@ -3,8 +3,8 @@
\textbf{Evaluation criteria:}
Each project has three elements that are graded: (i) the code,
(ii) the quality of the figures, and (iii) the presentation (see below).
You can view the evaluation criteria on the
\emph{SciCompScoreSheet.pdf} on Ilias.
\vspace{1ex}
\textbf{Dates:}
@ -18,38 +18,38 @@
\vspace{1ex}
\textbf{Files:}
Please hand in your presentation as a pdf file. Bundle
Hand in your presentation as a pdf file. Bundle
everything (the pdf, the code, and the data) into a {\em single}
zip-file.
Hint: make the zip file you want to upload, unpack it somewhere
else and check if your main script is still running properly.
\vspace{1ex}
\textbf{Code:}
The code should be executable without any further
adjustments from our side. A single \texttt{main.m} script
should coordinate the analysis by calling functions and
sub-scripts and should produce the {\em same} figures
(\texttt{saveas()}-function, pdf or png format) that you use in
your slides. The code should be properly commented and
The code must be executable without any further adjustments from
our side. (Test it!) A single \texttt{main.m} script coordinates
the analysis by calling functions and sub-scripts and produces
the {\em same} figures (\texttt{saveas()}-function, pdf or png
format) that you use in your slides. The code must be
comprehensible by a third person (use proper and consistent
variable and function names).
% Hint: make the zip file you want to upload, unpack it
% somewhere else and check if your main script is running.
\vspace{1ex}
\emph{Please write your name and matriculation number as a
\emph{Please note your name and matriculation number as a
comment at the top of the \texttt{main.m} script.}
\vspace{1ex}
\textbf{Presentation:}
The presentation should be {\em at most} 10min long and be held
in English. In the presentation you should present figures
introducing, explaining, showing, and discussing your data,
methods, and results. All data-related figures you show in the
presentation should be produced by your program --- no editing
presentation must be produced by your program --- no editing
or labeling by PowerPoint or other software. It is always a good
idea to illustrate the problem with basic plots of the
raw-data. Make sure the axis labels are large enough!

11
projects/project.mk
View File

@ -1,14 +1,10 @@
BASENAME=$(subst project_,,$(notdir $(CURDIR)))
latex:
pdflatex $(BASENAME).tex
pdflatex $(BASENAME).tex
pdf: $(BASENAME).pdf
$(BASENAME).pdf : $(BASENAME).tex ../header.tex ../instructions.tex
pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true
pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get" && pdflatex -interaction=scrollmode $< || true
watch :
@ -19,6 +15,11 @@ clean:
rm -rf *.log *.aux *.out auto
rm -f `basename *.tex .tex`.pdf
rm -f *.zip
pdflatex $(BASENAME).tex
latex:
pdflatex $(BASENAME).tex
pdflatex $(BASENAME).tex
zip: latex
rm -f zip $(BASENAME).zip

45
projects/project_eod/eod.tex
View File

@ -11,31 +11,48 @@
%%%%%%%%%%%%%% Questions %%%%%%%%%%%%%%%%%%%%%%%%%
Weakly electric fish employ their self-generated electric field for
prey-capture, navigation and also communication. In many of these fish
the {\em electric organ discharge} (EOD) is well described by a
combination of a sine-wave and a few of its harmonics (integer
multiples of the fundamental frequency).
\begin{questions}
\question In the data file {\tt EOD\_data.mat} you find a time trace
and the {\em electric organ discharge (EOD)} of a weakly electric
fish {\em Apteronotus leptorhynchus}.
\question In the data file {\tt EOD\_data.mat} you find two
variables. The first contains the time at which the EOD was sampled
and the second the acutal EOD recording of a weakly electric fisch
of the species {\em Apteronotus leptorhynchus}.
\begin{parts}
\part Load and plot the data in an appropriate way. Time is in
seconds and the voltage is in mV/cm.
\part Fit the following curve to the eod (select a small time
window, containing only 2 or three electric organ discharges, for
\part Load the data and create a plot showing the data. Time is given in
seconds and the voltage is given in mV/cm.
\part Fit the following curve to the EOD (select a \textbf{small} time
window, containing only two or three electric organ discharges, for
fitting, not the entire trace) using least squares:
$$f_{\omega_0,b_0,\varphi_1, ...,\varphi_n}(t) = b_0 +
\sum_{j=1}^n \alpha_j \cdot \sin(2\pi j\omega_0\cdot t + \varphi_j
).$$ $\omega_0$ is called {\em fundamental frequency}. The single
).$$ $\omega_0$ is called the {\em fundamental frequency}. The single
terms $\alpha_j \cdot \sin(2\pi j\omega_0\cdot t + \varphi_j )$
are called {\em harmonic components}. The variables $\varphi_j$
are called the {\em harmonic components}. The variables $\varphi_j$
are called {\em phases}, the $\alpha_j$ are the amplitudes. For
the beginning choose $n=3$.
\part Try different choices of $n$ and see how the fit
changes. Plot the fits and the original curve for different
choices of $n$. Also plot the fitting error as a function of
$n$.
changes. Plot the fits and the section of the original curve that
you used for fitting for different choices of $n$.
\part \label{fiterror} Plot the fitting error as a function of $n$.
What do you observe?
\part Another way to quantify the quality of the fit is to compute
the correlation coefficient between the fit and the
data. Illustrate this correlation for a few values of $n$. Plot
the correlation coefficient as a function of $n$. What is the
minimum $n$ needed for a good fit? How does this compare to the
results from (\ref{fiterror})?
\part Plot the amplitudes $\alpha_j$ and phases $\varphi_j$ as a
function of the frequencies $\omega_j$ --- the amplitude and phase
spectra, also called ``Bode plot''.
\part Why does the fitting fail when you try to fit the entire recording?
\part (optional) If you want you can also play the different fits
and the original as sound.
and the original as sound (check the help).
\end{parts}
\end{questions}

2
projects/project_eyetracker/eyetracker.tex
View File

@ -33,7 +33,7 @@ shifts. The eye movements during training and test are recorded.
speed and/or accelerations.
\part Detect and correct the eye traces for instances in which the
eye was not correctly detected. Interpolate linearily in these sections.
\part Create a 'heatmap' plot that shows the eye trajectories
\part Create a 'heatmap' plot of the eye-positions
for one or two (nice) trials.
\part Use the \verb+kmeans+ clustering function to
identify fixation points. Manually select a good number of cluster

16
projects/project_fano_slope/fano_slope.tex
View File

@ -1,6 +1,6 @@
\documentclass[a4paper,12pt,pdftex]{exam}
\newcommand{\ptitle}{Stimulus discrimination}
\newcommand{\ptitle}{Stimulus discrimination: gain}
\input{../header.tex}
\firstpagefooter{Supervisor: Jan Benda}{phone: 29 74573}%
{email: jan.benda@uni-tuebingen.de}
@ -95,10 +95,11 @@ spikes = lifboltzmanspikes(trials, input, tmax, gain);
the neuron? Plot them for the four different values of the gain
used in (a).
\part Think about a measure based on the spike-count histograms
that quantifies how well the two stimuli can be distinguished
based on the spike counts. Plot the dependence of this measure as
a function of the gain of the neuron.
\part \label{discrmeasure} Think about a measure based on the
spike-count histograms that quantifies how well the two stimuli
can be distinguished based on the spike counts. Plot the
dependence of this measure as a function of the gain of the
neuron.
%
For which gains can the two stimuli perfectly discriminated?
@ -110,6 +111,11 @@ spikes = lifboltzmanspikes(trials, input, tmax, gain);
results in the best discrimination performance. How can you
quantify ``best discrimination'' performance?
\part Another way to quantify the discriminability of the spike
counts in response to the two stimuli is to apply an appropriate
statistical test and check for significant differences. How does
this compare to your findings from (\ref{discrmeasure})?
\end{parts}
\end{questions}

16
projects/project_fano_time/fano_time.tex
View File

@ -1,6 +1,6 @@
\documentclass[a4paper,12pt,pdftex]{exam}
\newcommand{\ptitle}{Stimulus discrimination}
\newcommand{\ptitle}{Stimulus discrimination: time}
\input{../header.tex}
\firstpagefooter{Supervisor: Jan Benda}{phone: 29 74573}%
{email: jan.benda@uni-tuebingen.de}
@ -87,10 +87,11 @@ input = 15.0; % I_2
observation time $T$? Plot them for four different values of $T$
(use values of 10\,ms, 100\,ms, 300\,ms and 1\,s).
\part Think about a measure based on the spike-count histograms
that quantifies how well the two stimuli can be distinguished
based on the spike counts. Plot the dependence of this measure as
a function of the observation time $T$.
\part \label{discrmeasure} Think about a measure based on the
spike-count histograms that quantifies how well the two stimuli
can be distinguished based on the spike counts. Plot the
dependence of this measure as a function of the observation time
$T$.
For which observation times can the two stimuli perfectly
discriminated?
@ -103,6 +104,11 @@ input = 15.0; % I_2
results in the best discrimination performance. How can you
quantify ``best discrimination'' performance?
\part Another way to quantify the discriminability of the spike
counts in response to the two stimuli is to apply an appropriate
statistical test and check for significant differences. How does
this compare to your findings from (\ref{discrmeasure})?
\end{parts}
\end{questions}

71
projects/project_isicorrelations/isicorrelations.tex
View File

@ -1,6 +1,6 @@
\documentclass[a4paper,12pt,pdftex]{exam}
\newcommand{\ptitle}{Interspike-intervall correlations}
\newcommand{\ptitle}{Adaptation and interspike-interval correlations}
\input{../header.tex}
\firstpagefooter{Supervisor: Jan Benda}{phone: 29 74573}%
{email: jan.benda@uni-tuebingen.de}
@ -20,10 +20,10 @@
Explore the dependence of interspike interval correlations on the firing rate,
adaptation time constant and noise level.
The neuron is a neuron with an adaptation current.
It is implemented in the file \texttt{lifadaptspikes.m}. Call it
with the following parameters:
\begin{lstlisting}
The neuron is a neuron with an adaptation current. It is
implemented in the file \texttt{lifadaptspikes.m}. Call it with the
following parameters:
\begin{lstlisting}
trials = 10;
tmax = 50.0;
input = 10.0; % the input I
@ -31,7 +31,7 @@ Dnoise = 1e-2; % noise strength
adapttau = 0.1; % adaptation time constant in seconds
adaptincr = 0.5; % adaptation strength
spikes = lifadaptspikes( trials, input, tmax, Dnoise, adapttau, adaptincr );
spikes = lifadaptspikes(trials, input, tmax, Dnoise, adapttau, adaptincr);
\end{lstlisting}
The returned \texttt{spikes} is a cell array with \texttt{trials} elements, each being a vector
of spike times (in seconds) computed for a duration of \texttt{tmax} seconds.
@ -39,29 +39,42 @@ spikes = lifadaptspikes( trials, input, tmax, Dnoise, adapttau, adaptincr );
and the adaptation time constant via \texttt{adapttau}.
\begin{parts}
\part Measure the intensity-response curve of the neuron, that is the mean firing rate
as a function of the input for a range of inputs from 0 to 120.
\part Compute the correlations between each interspike interval $T_i$ and the next one $T_{i+1}$
(serial interspike interval correlation at lag 1). Plot this correlation as a function of the
firing rate by varying the input as in (a).
\part How does this dependence change for different values of the adaptation
time constant \texttt{adapttau}? Use values between 10\,ms and
1\,s for \texttt{adapttau}.
\part Determine the firing rate at which the minimum interspike interval correlation
occurs. How does the minimum correlation and this firing rate
depend on the adaptation time constant \texttt{adapttau}?
\part How do the results change if the level of the intrinsic noise \texttt{Dnoise} is modified?
Use values of 1e-4, 1e-3, 1e-2, 1e-1, and 1 for \texttt{Dnoise}.
\uplevel{If you still have time you can continue with the following question:}
\part How do the interspike interval distributions look like for the different noise levels
at some example values for the input and the adaptation time constant?
\part Show a spike-raster plot and a time-resolved firing rate of
the neuron for an input current of 50 for three different
adaptation time constants \texttt{adapttau} (10\,m, 100\,ms,
1\,s). How do the neural responses look like and how do they
depend on the adaptation time constant?
\uplevel{For all the following analysis we only use the spike
times of the steady-state response, i.e. we skip all spikes
occuring before at least three times the adaptation time
constant.}
\part \label{ficurve} Measure the intensity-response curve of the
neuron, that is the mean firing rate as a function of the input
for a range of inputs from 0 to 120.
\part Additionally compute the correlations between each
interspike interval $T_i$ and the next one $T_{i+1}$ (serial
interspike interval correlation at lag 1) for the same range of
inputs as in (\ref{ficurve}). Plot the correlation as a function
of the input.
\part How does the intensity-response curve and the
interspike-interval correlations depend on the adaptation time
constant \texttt{adapttau}? Use several values between 10\,ms and
1\,s for \texttt{adapttau} (logarithmically distributed).
\part Determine the firing rate at which the minimum interspike
interval correlation occurs. How does the minimum correlation and
this firing rate (or the inverse of it, the mean interspike
interval) depend on the adaptation time constant
\texttt{adapttau}? Is this dependence siginificant? If yes, can
you explain this dependence?
\part How do all the results change if the level of the intrinsic
noise \texttt{Dnoise} is modified? Use values of 1e-4, 1e-3,
1e-2, 1e-1, and 1 for \texttt{Dnoise}.
\end{parts}

7
projects/project_isicorrelations/lifadaptspikes.m
View File

@ -41,11 +41,8 @@ function spikes = lifadaptspikes( trials, input, tmaxdt, D, tauadapt, adaptincr
if v >= vthresh
v = vreset;
a = a + adaptincr/tauadapt;
spiketime = i*dt;
if spiketime > 4.0*tauadapt
times(j) = spiketime - 4.0*tauadapt;
j = j + 1;
end
times(j) = i*dt;
j = j + 1;
end
end
spikes{k} = times;

11
projects/project_isicorrelations/solution/firingrate.m Normal file
View File

@ -0,0 +1,11 @@
function [ rate ] = firingrate(spikes, t0, t1)
% compute the firing rate from spikes between t0 and t1
rates = zeros(length(spikes), 1);
for k = 1:length(spikes)
spiketimes = spikes{k};
rates(k) = length(spiketimes((spiketimes>=t0)&(spiketimes<=t1)))/(t1-t0);
end
rate = mean(rates);
end

36
projects/project_isicorrelations/solution/isicorrelations.m Normal file
View File

@ -0,0 +1,36 @@
trials = 5;
tmax = 10.0;
Dnoise = 1e-2; % noise strength
adapttau = 0.1; % adaptation time constant in seconds
adaptincr = 0.5; % adaptation strength
t0=2.0;
t1=tmax;
maxlag = 5;
taus = [0.01, 0.1, 1.0];
colors = ['r', 'b', 'g'];
for j = 1:length(taus)
adapttau = taus(j);
% f-I curves:
Is = [0:10:80];
rate = zeros(length(Is),1);
corr = zeros(length(Is),maxlag);
for k = 1:length(Is)
input = Is(k);
spikes = lifadaptspikes(trials, input, tmax, Dnoise, adapttau, adaptincr);
rate(k) = firingrate(spikes, t0, t1);
corr(k,:) = serialcorr(spikes, t0, t1, maxlag);
end
figure(1);
hold on
plot(Is, rate, colors(j));
hold off
figure(2);
hold on
plot(Is, corr(:,2), colors(j));
hold off
end
pause

50
projects/project_isicorrelations/solution/lifadaptspikes.m Normal file
View File

@ -0,0 +1,50 @@
function spikes = lifadaptspikes( trials, input, tmaxdt, D, tauadapt, adaptincr )
% Generate spike times of a leaky integrate-and-fire neuron
% with an adaptation current
% trials: the number of trials to be generated
% input: the stimulus either as a single value or as a vector
% tmaxdt: in case of a single value stimulus the duration of a trial
% in case of a vector as a stimulus the time step
% D: the strength of additive white noise
% tauadapt: adaptation time constant
% adaptincr: adaptation strength
tau = 0.01;
if nargin < 4
D = 1e0;
end
if nargin < 5
tauadapt = 0.1;
end
if nargin < 6
adaptincr = 1.0;
end
vreset = 0.0;
vthresh = 10.0;
dt = 1e-4;
if max( size( input ) ) == 1
input = input * ones( ceil( tmaxdt/dt ), 1 );
else
dt = tmaxdt;
end
spikes = cell( trials, 1 );
for k=1:trials
times = [];
j = 1;
v = vreset + (vthresh-vreset)*rand();
a = 0.0;
noise = sqrt(2.0*D)*randn( length( input ), 1 )/sqrt(dt);
for i=1:length( noise )
v = v + ( - v - a + noise(i) + input(i))*dt/tau;
a = a + ( - a )*dt/tauadapt;
if v >= vthresh
v = vreset;
a = a + adaptincr/tauadapt;
times(j) = i*dt;
j = j + 1;
end
end
spikes{k} = times;
end
end

15
projects/project_isicorrelations/solution/serialcorr.m Normal file
View File

@ -0,0 +1,15 @@
function [ isicorrs ] = serialcorr(spikes, t0, t1, maxlag)
% compute the serial correlations from spikes between t0 and t1
ics = zeros(length(spikes), maxlag);
for k = 1:length(spikes)
spiketimes = spikes{k};
isis = diff(spiketimes((spiketimes>=t0)&(spiketimes<=t1)));
if length(isis) > 2*maxlag
for j=1:maxlag
ics(k, j) = corr(isis(j:end)', isis(1:end+1-j)');
end
end
end
isicorrs = mean(ics, 1);
end

6
projects/project_lif/lif.tex
View File

@ -96,7 +96,7 @@ time = [0.0:dt:tmax]; % t_i
Write a function that implements this leaky integrate-and-fire
neuron by expanding the function for the passive neuron
appropriate. The function returns a vector of spike times.
appropriately. The function returns a vector of spike times.
Illustrate how this model works by appropriate plot(s) and
input(s) $E(t)$, e.g. constant inputs lower and higher than the
@ -115,8 +115,8 @@ time = [0.0:dt:tmax]; % t_i
r = \frac{n-1}{t_n - t_1}
\end{equation}
What do you observe? Does the firing rate encode the frequency of
the stimulus? Look at the spike trains in response the sine waves
to figure out what is going on.
the stimulus? Look at the spike trains in response to the sine
waves to figure out what is going on.
\end{parts}
\end{questions}

6
projects/project_noiseficurves/noiseficurves.tex
View File

@ -63,8 +63,8 @@ spikes = lifspikes(trials, current, tmax, Dnoise);
of this neuron?
\part Compute the $f$-$I$ curves of neurons with various noise
strengths \texttt{Dnoise}. Use for example $D_{noise} = 1e-3$,
$1e-2$, and $1e-1$.
strengths \texttt{Dnoise}. Use for example $D_{noise} = 10^{-3}$,
$10^{-2}$, and $10^{-1}$.
How does the intrinsic noise influence the response curve?
@ -76,6 +76,8 @@ spikes = lifspikes(trials, current, tmax, Dnoise);
responses of the four different neurons to the same input, or by
the same resulting mean firing rate.
How do the responses differ?
\part Let's now use as an input to the neuron a 1\,s long sine
wave $I(t) = I_0 + A \sin(2\pi f t)$ with offset current $I_0$,
amplitude $A$, and frequency $f$. Set $I_0=5$, $A=4$, and

5
projects/project_pca_natural_images/pca_natural_images.tex
View File

@ -32,9 +32,8 @@ In you zip file you find a natural image called {\tt natimg.jpg}.
\begin{thebibliography}{1}
\bibitem{BG} Buchsbaum, G., \& Gottschalk, A. (1983). Trichromacy,
opponent colours coding and optimum colour information transmission
in the retina. Proceedings of the Royal Society of London. Series B,
Containing Papers of a Biological Character. Royal Society (Great
Britain), 220(1218), 89–113.
in the retina. Proceedings of the Royal Society of London B. Royal
Society (Great Britain), 220(1218), 89–113.
\end{thebibliography}

50
projects/project_photoreceptor/photoreceptor.tex
View File

@ -11,11 +11,16 @@
%%%%%%%%%%%%%% Questions %%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Analysis of insect photoreceptor data.}
\section*{Light responses of an insect photoreceptor.}
In this project you will analyse data from intracellular recordings of
a fly R\.1--6 photoreceptor. The membrane potential of the
photoreceptor was recorded while the cell was stimulated with a
light stimulus.
a fly R\,1--6 photoreceptor. These cells show graded membrane
potential changes in response to a light stimulus. The membrane
potential of the photoreceptor was recorded while the cell was
stimulated with a light stimulus. Intracellular recordings often
suffer from drifts in the resting potential. This leads to a large
variability in the responses which is technical and not a cellular
property. To compensate for such drifts trials are aligned to the
resting potential before stimulus onset.
\begin{questions}
\question{} The accompanying dataset (photoreceptor\_data.zip)
@ -24,23 +29,36 @@ light stimulus.
\textit{voltage} a matrix with the recorded membrane potential from
10 consecutive trials, (ii) \textit{time} a matrix with the
time-axis for each trial, and (iii) \textit{trace\_meta} a structure
that stores several metadata. This is the place where you find the
\emph{amplitude}, that is the voltage that drives the light
stimulus, i.e. the light-intensity.
that stores several metadata including the \emph{amplitude} value
that is the voltage used to drive the light stimulus. (Note that
this voltage is only a proxy for the true light intensity. Twice the
voltage does not lead to twice the light intensity. Within this
project, however, you can treat it as if it was the intensity.)
\begin{parts}
\part{} Create a plot of the raw data. Plot the average response as
a function of time. This plot should also show the
across-trial variability.\\[0.5ex]
\part{} You will notice that the responses have three main parts, a
\part Create a plot of the raw data. For each light intensity plot
the average response as a function of time. This plot should also
depict the across-trial variability in an appropriate way.
\part You will notice that the responses have three main parts, a
pre-stimulus phase, the phase in which the light was on, and
finally a post-stimulus phase. Create an characteristic curve that
plots the response strength as a function of the stimulus
intensity for the ``onset'' and the ``steady state''
phases.\\[0.5ex]
\part{} The light switches on at time zero. Estimate the delay between stimulus.\\[0.5ex]
\part{} You may also decide to analyze the post-stimulus response in some
more detail.
phases of the light response.
\part The light switches on at time zero. Estimate the delay
between stimulus and response.
\part Analyze the across trial variability in the ``onset'' and
``steady state''. Check for statistically significant differences.
\part The membrane potential shows some fluctuations (noise)
compare the noise before stimulus onset and in the steady state
phase of the response.
\part (optional) You may also analyze the post-stimulus response
in some more detail.
\end{parts}
\end{questions}

32
projects/project_stimulus_reconstruction/stimulus_reconstruction.tex
View File

@ -10,29 +10,29 @@
\input{../instructions.tex}
%%%%%%%%%%%%%% Questions %%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Reverse reconstruction of the stimulus evoking neuronal responses.}
\section*{Reverse reconstruction of the stimulus that evoked a neuronal response.}
To analyse encoding properties of a neuron one often calculates the
Spike-Triggered-Average (STA).
\[ STA(\tau) = \frac{1}{\langle n \rangle} \left\langle
\displaystyle\sum_{i=1}^{n}{s(t_i - \tau)} \right\rangle \]
The STA is the average stimulus that led to a spike in the neuron and
is calculated by cutting out snippets form the stimulus centered on
the respective spike time. The Spike-Triggered-Average can be used to
reconstruct the stimulus a neuron has been stimulated with.
Spike-Triggered-Average (STA). The STA is the average stimulus that
led to a spike in the neuron and is calculated by cutting out snippets
form the stimulus centered on the respective spike time:
\[ STA(\tau) = \frac{1}{n} \displaystyle\sum_{i=1}^{n}{s(t_i - \tau)} \],
where $n$ is the number of trials and $t_i$ is the time of the
$i_{th}$ spike. The Spike-Triggered-Average can be used to reconstruct
the stimulus from the neuronal response. The reconstructed stimulus
can then be compared to the original stimulus.
\begin{questions}
\question In the accompanying files you find the spike responses of
P-units and pyramidal neurons of the weakly electric fish
a p-type electroreceptor afferent (P-unit) and a pyramidal neurons
recorded in the hindbrain of the weakly electric fish
\textit{Apteronotus leptorhynchus}. The respective stimuli are
stored in separate files. The data is sampled with 20\,kHz temporal
resolution and spike times are given in seconds. Start with the
P-unit and, in the end, apply the same functions to the pyramidal
data.
P-unit and, in the end, apply the same analyzes/functions to the
pyramidal data.
\begin{parts}
\part Estimate the STA and plot it.
\part Implement a function that does the reconstruction of the
stimulus using the STA.
\part Implement a function that does the reverse reconstruction and uses the STA to recopnstruct the stimulus.
\part Implement a function that estimates the reconstruction
error using the mean-square-error and express it relative to the
variance of the original stimulus.
@ -44,8 +44,8 @@ reconstruct the stimulus a neuron has been stimulated with.
\part Analyze the robustness of the reconstruction: Estimate
the STA with less and less data and estimate the reconstruction
error.
\part Plot the reconstruction error as a function of the data
amount used to estimate the STA.
\part Plot the reconstruction error as a function of the amount of data
used to estimate the STA.
\part Repeat the above steps for the pyramidal neuron, what do you
observe?
\end{parts}

34
regression/exercises/exercises01.tex
View File

@ -20,7 +20,7 @@
\else
\newcommand{\stitle}{}
\fi
\header{{\bfseries\large Exercise 11\stitle}}{{\bfseries\large Gradient descent}}{{\bfseries\large January 9th, 2018}}
\header{{\bfseries\large Exercise 10\stitle}}{{\bfseries\large Gradient descent}}{{\bfseries\large December 17th, 2018}}
\firstpagefooter{Dr. Jan Grewe}{Phone: 29 74588}{Email:
jan.grewe@uni-tuebingen.de}
\runningfooter{}{\thepage}{}
@ -58,22 +58,24 @@
\begin{questions}
\question Implement the gradient descent for finding the parameters
of a straigth line that we want to fit to the data in the file
\emph{lin\_regression.mat}.
\question We want to fit the straigth line \[ y = mx+b \] to the
data in the file \emph{lin\_regression.mat}.
In the lecture we already prepared most of the necessary functions:
1. the error function (\code{meanSquareError()}), 2. the cost
function (\code{lsqError()}), and 3. the gradient
(\code{lsqGradient()}). Read chapter 8 ``Optimization and gradient
descent'' in the script, in particular section 8.4 and exercise 8.4!
In the lecture we already prepared the cost function
(\code{lsqError()}), and the gradient (\code{lsqGradient()}) (read
chapter 8 ``Optimization and gradient descent'' in the script, in
particular section 8.4 and exercise 8.4!). With these functions in
place we here want to implement a gradient descend algorithm that
finds the minimum of the cost function and thus the slope and
intercept of the straigth line that minimizes the squared distance
to the data values.
The algorithm for the descent towards the minimum of the cost
function is as follows:
\begin{enumerate}
\item Start with some arbitrary parameter values $\vec p_0 = (m_0, b_0)$
for the slope and the intercept of the straight line.
\item Start with some arbitrary parameter values (intercept $b_0$
and slope $m_0$, $\vec p_0 = (b_0, m_0)$ for the slope and the
intercept of the straight line.
\item \label{computegradient} Compute the gradient of the cost function
at the current values of the parameters $\vec p_i$.
\item If the magnitude (length) of the gradient is smaller than some
@ -106,9 +108,11 @@
\lstinputlisting{../code/descentfit.m}
\end{solution}
\part Find the position of the minimum of the cost function by
means of the \code{min()} function. Compare with the result of the
gradient descent method. Vary the value of $\epsilon$ and the
\part For checking the gradient descend method from (a) compare
its result for slope and intercept with the position of the
minimum of the cost function that you get when computing the cost
function for many values of the slope and intercept and then using
the \code{min()} function. Vary the value of $\epsilon$ and the
minimum gradient. What are good values such that the gradient
descent gets closest to the true minimum of the cost function?
\begin{solution}

4
regression/lecture/regression-chapter.tex
View File

@ -16,6 +16,10 @@
\input{regression}
\section{Improvements}
Adapt function arguments to matlabs polyfit. That is: first the data
(x,y) and then the parameter vector. p(1) is slope, p(2) is intercept.
\section{Fitting in practice}
Fit with matlab functions lsqcurvefit, polyfit