diff --git a/likelihood/exercises/exercises01-de.tex b/likelihood/exercises/exercises01-de.tex new file mode 100644 index 0000000..1c3dcdf --- /dev/null +++ b/likelihood/exercises/exercises01-de.tex @@ -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} \ No newline at end of file diff --git a/likelihood/exercises/exercises01.tex b/likelihood/exercises/exercises01.tex index 1c3dcdf..ca715f5 100644 --- a/likelihood/exercises/exercises01.tex +++ b/likelihood/exercises/exercises01.tex @@ -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} diff --git a/likelihood/exercises/instructions.tex b/likelihood/exercises/instructions.tex index 12e9e24..7fe599d 100644 --- a/likelihood/exercises/instructions.tex +++ b/likelihood/exercises/instructions.tex @@ -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 diff --git a/likelihood/exercises/mlepropest.m b/likelihood/exercises/mlepropest.m new file mode 100644 index 0000000..352c28c --- /dev/null +++ b/likelihood/exercises/mlepropest.m @@ -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); diff --git a/likelihood/exercises/mlepropest.pdf b/likelihood/exercises/mlepropest.pdf new file mode 100644 index 0000000..0a22974 Binary files /dev/null and b/likelihood/exercises/mlepropest.pdf differ diff --git a/likelihood/exercises/mlestd.m b/likelihood/exercises/mlestd.m index de37a56..f28d39e 100644 --- a/likelihood/exercises/mlestd.m +++ b/likelihood/exercises/mlestd.m @@ -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); diff --git a/likelihood/exercises/mlestd.pdf b/likelihood/exercises/mlestd.pdf index f01f927..5246eda 100644 Binary files a/likelihood/exercises/mlestd.pdf and b/likelihood/exercises/mlestd.pdf differ diff --git a/linearalgebra/code/coordinaterafo2.m b/linearalgebra/code/coordinaterafo2.m deleted file mode 100644 index 054ac6d..0000000 --- a/linearalgebra/code/coordinaterafo2.m +++ /dev/null @@ -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 diff --git a/linearalgebra/code/coordinatetrafo.m b/linearalgebra/code/coordinatetrafo.m index 7c8507c..de62642 100644 --- a/linearalgebra/code/coordinatetrafo.m +++ b/linearalgebra/code/coordinatetrafo.m @@ -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; diff --git a/linearalgebra/code/covareigen.m b/linearalgebra/code/covareigen.m index cd6025f..03431f4 100644 --- a/linearalgebra/code/covareigen.m +++ b/linearalgebra/code/covareigen.m @@ -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' ); diff --git a/linearalgebra/code/covareigen3.m b/linearalgebra/code/covareigen3.m index 636d9cb..ec4bee4 100644 --- a/linearalgebra/code/covareigen3.m +++ b/linearalgebra/code/covareigen3.m @@ -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' ); diff --git a/linearalgebra/code/covareigen3examples.m b/linearalgebra/code/covareigen3examples.m index f27c76c..faefbde 100644 --- a/linearalgebra/code/covareigen3examples.m +++ b/linearalgebra/code/covareigen3examples.m @@ -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; diff --git a/linearalgebra/code/covareigenexamples.m b/linearalgebra/code/covareigenexamples.m index 2f7a480..5b7262d 100644 --- a/linearalgebra/code/covareigenexamples.m +++ b/linearalgebra/code/covareigenexamples.m @@ -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 ); diff --git a/linearalgebra/code/covarmatrix.m b/linearalgebra/code/covarmatrix.m new file mode 100644 index 0000000..ff4368e --- /dev/null +++ b/linearalgebra/code/covarmatrix.m @@ -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() diff --git a/linearalgebra/code/covarmatrixeigen.m b/linearalgebra/code/covarmatrixeigen.m new file mode 100644 index 0000000..c0a6f74 --- /dev/null +++ b/linearalgebra/code/covarmatrixeigen.m @@ -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() diff --git a/linearalgebra/code/matrixbox.m b/linearalgebra/code/matrixbox.m index 6cc0aa1..4686dba 100644 --- a/linearalgebra/code/matrixbox.m +++ b/linearalgebra/code/matrixbox.m @@ -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 diff --git a/linearalgebra/code/pca2d.m b/linearalgebra/code/pca2d.m index 6879f8a..64c3d21 100644 --- a/linearalgebra/code/pca2d.m +++ b/linearalgebra/code/pca2d.m @@ -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' ); diff --git a/linearalgebra/code/pca2dexamples.m b/linearalgebra/code/pca2dexamples.m index 252ceda..52ae0b9 100644 --- a/linearalgebra/code/pca2dexamples.m +++ b/linearalgebra/code/pca2dexamples.m @@ -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 diff --git a/linearalgebra/code/pca3d.m b/linearalgebra/code/pca3d.m index 6f2a6bd..c485638 100644 --- a/linearalgebra/code/pca3d.m +++ b/linearalgebra/code/pca3d.m @@ -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 + diff --git a/linearalgebra/code/pca3dexamples.m b/linearalgebra/code/pca3dexamples.m index d2d790f..5637c1a 100644 --- a/linearalgebra/code/pca3dexamples.m +++ b/linearalgebra/code/pca3dexamples.m @@ -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 diff --git a/linearalgebra/code/simplematrixbox.m b/linearalgebra/code/simplematrixbox.m index bb94254..b96f5cd 100644 --- a/linearalgebra/code/simplematrixbox.m +++ b/linearalgebra/code/simplematrixbox.m @@ -17,5 +17,5 @@ function simplematrixbox( a, s ) xlim( [-2 2 ] ) ylim( [-2 2] ) title( s ); - waitforbuttonpress; + pause(1.0); end diff --git a/linearalgebra/code/spikesortingwave.m b/linearalgebra/code/spikesortingwave.m index da7c251..ce2ebec 100644 --- a/linearalgebra/code/spikesortingwave.m +++ b/linearalgebra/code/spikesortingwave.m @@ -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 $(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) diff --git a/linearalgebra/exercises/UT_WBMW_Black_RGB.pdf b/linearalgebra/exercises/UT_WBMW_Black_RGB.pdf new file mode 100644 index 0000000..9aed921 Binary files /dev/null and b/linearalgebra/exercises/UT_WBMW_Black_RGB.pdf differ diff --git a/linearalgebra/exercises/covariance.m b/linearalgebra/exercises/covariance.m new file mode 100644 index 0000000..2f60060 --- /dev/null +++ b/linearalgebra/exercises/covariance.m @@ -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); diff --git a/linearalgebra/exercises/covariance.pdf b/linearalgebra/exercises/covariance.pdf new file mode 100644 index 0000000..6eb41e8 Binary files /dev/null and b/linearalgebra/exercises/covariance.pdf differ diff --git a/linearalgebra/exercises/exercises01.tex b/linearalgebra/exercises/exercises01.tex new file mode 100644 index 0000000..4264a75 --- /dev/null +++ b/linearalgebra/exercises/exercises01.tex @@ -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} \ No newline at end of file diff --git a/linearalgebra/exercises/extdata.mat b/linearalgebra/exercises/extdata.mat new file mode 100644 index 0000000..619c0c1 Binary files /dev/null and b/linearalgebra/exercises/extdata.mat differ diff --git a/linearalgebra/exercises/instructions.tex b/linearalgebra/exercises/instructions.tex new file mode 100644 index 0000000..3041d3e --- /dev/null +++ b/linearalgebra/exercises/instructions.tex @@ -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} diff --git a/linearalgebra/exercises/nonlinearcorrelation.m b/linearalgebra/exercises/nonlinearcorrelation.m new file mode 100644 index 0000000..bcf9a59 --- /dev/null +++ b/linearalgebra/exercises/nonlinearcorrelation.m @@ -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); diff --git a/linearalgebra/exercises/nonlinearcorrelation.pdf b/linearalgebra/exercises/nonlinearcorrelation.pdf new file mode 100644 index 0000000..b65389b Binary files /dev/null and b/linearalgebra/exercises/nonlinearcorrelation.pdf differ diff --git a/linearalgebra/exercises/pca2d-2.pdf b/linearalgebra/exercises/pca2d-2.pdf new file mode 100644 index 0000000..f413a21 Binary files /dev/null and b/linearalgebra/exercises/pca2d-2.pdf differ diff --git a/linearalgebra/exercises/pca2d.m b/linearalgebra/exercises/pca2d.m new file mode 100644 index 0000000..972590b --- /dev/null +++ b/linearalgebra/exercises/pca2d.m @@ -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 diff --git a/linearalgebra/exercises/pca3d-2.pdf b/linearalgebra/exercises/pca3d-2.pdf new file mode 100644 index 0000000..57b6225 Binary files /dev/null and b/linearalgebra/exercises/pca3d-2.pdf differ diff --git a/linearalgebra/exercises/pca3d.m b/linearalgebra/exercises/pca3d.m new file mode 100644 index 0000000..4bd05dc --- /dev/null +++ b/linearalgebra/exercises/pca3d.m @@ -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 diff --git a/linearalgebra/exercises/savefigpdf.m b/linearalgebra/exercises/savefigpdf.m new file mode 100644 index 0000000..2a11c11 --- /dev/null +++ b/linearalgebra/exercises/savefigpdf.m @@ -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 + diff --git a/linearalgebra/exercises/spikesorting.m b/linearalgebra/exercises/spikesorting.m new file mode 100644 index 0000000..af958ca --- /dev/null +++ b/linearalgebra/exercises/spikesorting.m @@ -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 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). \ No newline at end of file +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 \ No newline at end of file diff --git a/projects/evaluation.tex b/projects/evaluation.tex index 51a3071..91afda7 100644 --- a/projects/evaluation.tex +++ b/projects/evaluation.tex @@ -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 diff --git a/projects/instructions.tex b/projects/instructions.tex index d68dd46..d4791cb 100644 --- a/projects/instructions.tex +++ b/projects/instructions.tex @@ -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! diff --git a/projects/project.mk b/projects/project.mk index 2f8d3b9..1a9bf0b 100644 --- a/projects/project.mk +++ b/projects/project.mk @@ -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 diff --git a/projects/project_eod/eod.tex b/projects/project_eod/eod.tex index 647a756..59b7254 100644 --- a/projects/project_eod/eod.tex +++ b/projects/project_eod/eod.tex @@ -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} diff --git a/projects/project_eyetracker/eyetracker.tex b/projects/project_eyetracker/eyetracker.tex index 9ad9afb..017c368 100644 --- a/projects/project_eyetracker/eyetracker.tex +++ b/projects/project_eyetracker/eyetracker.tex @@ -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 diff --git a/projects/project_fano_slope/fano_slope.tex b/projects/project_fano_slope/fano_slope.tex index fde346d..a204e9a 100644 --- a/projects/project_fano_slope/fano_slope.tex +++ b/projects/project_fano_slope/fano_slope.tex @@ -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} diff --git a/projects/project_fano_time/fano_time.tex b/projects/project_fano_time/fano_time.tex index e2960af..42f788a 100644 --- a/projects/project_fano_time/fano_time.tex +++ b/projects/project_fano_time/fano_time.tex @@ -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} diff --git a/projects/project_isicorrelations/isicorrelations.tex b/projects/project_isicorrelations/isicorrelations.tex index f9eee80..7b49a02 100644 --- a/projects/project_isicorrelations/isicorrelations.tex +++ b/projects/project_isicorrelations/isicorrelations.tex @@ -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} diff --git a/projects/project_isicorrelations/lifadaptspikes.m b/projects/project_isicorrelations/lifadaptspikes.m index e2de8fc..2236a76 100644 --- a/projects/project_isicorrelations/lifadaptspikes.m +++ b/projects/project_isicorrelations/lifadaptspikes.m @@ -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; diff --git a/projects/project_isicorrelations/solution/firingrate.m b/projects/project_isicorrelations/solution/firingrate.m new file mode 100644 index 0000000..ea85254 --- /dev/null +++ b/projects/project_isicorrelations/solution/firingrate.m @@ -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 + diff --git a/projects/project_isicorrelations/solution/isicorrelations.m b/projects/project_isicorrelations/solution/isicorrelations.m new file mode 100644 index 0000000..49d7985 --- /dev/null +++ b/projects/project_isicorrelations/solution/isicorrelations.m @@ -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 diff --git a/projects/project_isicorrelations/solution/lifadaptspikes.m b/projects/project_isicorrelations/solution/lifadaptspikes.m new file mode 100644 index 0000000..2236a76 --- /dev/null +++ b/projects/project_isicorrelations/solution/lifadaptspikes.m @@ -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 diff --git a/projects/project_isicorrelations/solution/serialcorr.m b/projects/project_isicorrelations/solution/serialcorr.m new file mode 100644 index 0000000..589337c --- /dev/null +++ b/projects/project_isicorrelations/solution/serialcorr.m @@ -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 diff --git a/projects/project_lif/lif.tex b/projects/project_lif/lif.tex index 72a36ea..068b85e 100644 --- a/projects/project_lif/lif.tex +++ b/projects/project_lif/lif.tex @@ -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} diff --git a/projects/project_noiseficurves/noiseficurves.tex b/projects/project_noiseficurves/noiseficurves.tex index 3add8d8..98f0aad 100644 --- a/projects/project_noiseficurves/noiseficurves.tex +++ b/projects/project_noiseficurves/noiseficurves.tex @@ -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 diff --git a/projects/project_pca_natural_images/pca_natural_images.tex b/projects/project_pca_natural_images/pca_natural_images.tex index 0091bd7..c82c625 100644 --- a/projects/project_pca_natural_images/pca_natural_images.tex +++ b/projects/project_pca_natural_images/pca_natural_images.tex @@ -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} diff --git a/projects/project_photoreceptor/photoreceptor.tex b/projects/project_photoreceptor/photoreceptor.tex index bebfbb7..350341b 100644 --- a/projects/project_photoreceptor/photoreceptor.tex +++ b/projects/project_photoreceptor/photoreceptor.tex @@ -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} diff --git a/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex b/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex index 0730ab0..eae2513 100644 --- a/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex +++ b/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex @@ -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} diff --git a/regression/exercises/exercises01.tex b/regression/exercises/exercises01.tex index c6eadc1..a58e3da 100644 --- a/regression/exercises/exercises01.tex +++ b/regression/exercises/exercises01.tex @@ -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} diff --git a/regression/lecture/regression-chapter.tex b/regression/lecture/regression-chapter.tex index 39b8a29..7a76010 100644 --- a/regression/lecture/regression-chapter.tex +++ b/regression/lecture/regression-chapter.tex @@ -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