356 lines
16 KiB
TeX
356 lines
16 KiB
TeX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\chapter{Stochastic simulations}
|
|
\label{stochasticsimulationschapter}
|
|
\exercisechapter{Stochastic simulations}
|
|
|
|
The real power of computers for data analysis is the possibility to
|
|
run simulations. Experimental data of almost unlimited sample sizes
|
|
can be simulated in no time. This allows to explore basic concepts,
|
|
the ones we introduce in the following chapters and many more, with
|
|
well controlled data sets that are free of confounding pecularities of
|
|
real data sets. With simulated data we can also test our own analysis
|
|
functions. More importantly, by means of simulations we can explore
|
|
possible outcomes of our experiments before we even started the
|
|
experiment. We could even explore possible results for regimes that we
|
|
cannot test experimentally. How dynamical systems, like for example
|
|
predator-prey interactions or the activity of neurons or whole brains,
|
|
evolve in time is a central application for simulations. The advent of
|
|
computers at the second half of the twentieth century pushed the
|
|
exciting field of nonlinear dynamical systems forward. Conceptually,
|
|
many kinds of simulations are very simple and are implemented in a few
|
|
lines of code.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Random numbers}
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics{randomnumbers}
|
|
\titlecaption{\label{randomnumbersfig} Random numbers.} {Numerical
|
|
random number generators return sequences of numbers that are as
|
|
random as possible. The returned values approximate a given
|
|
probability distribution. Here, for example, a uniform
|
|
distribution between zero and one (top left). The serial
|
|
correlation (auto-correlation) of the returned sequences is zero
|
|
at any lag except for zero (bottom left). However, by setting a
|
|
seed, defining the state, or spawning the random number generator,
|
|
the exact same sequence can be generated (right).}
|
|
\end{figure}
|
|
|
|
At the heart of many simulations are random numbers that we get from
|
|
\enterm[random number generator]{random number generators}. These are
|
|
numerical algorithms that return sequences of numbers that appear to
|
|
be as random as possible (\figref{randomnumbersfig}). If we draw
|
|
random numbers using, for example, the \code{rand()} function, then
|
|
these numbers are indeed uniformly distributed and have a mean of one
|
|
half. Subsequent numbers are independent of each other, i.e. the
|
|
autocorrelation function is zero everywhere except at lag
|
|
zero. However, numerical random number generators have a period, after
|
|
which they repeat the exact same sequence of numbers. This
|
|
differentiates them from truely random numbers and hence they are
|
|
called \enterm[random number generator!pseudo]{pseudo random number
|
|
generators}. In rare cases this periodicity can induce problems in
|
|
simulations whenever more random numbers than the period of the random
|
|
number generator are used. Luckily, nowadays the periods of random
|
|
nunmber generators are very large. Periods are at least in the range
|
|
of $2^{64}$, that is about $10^{18}$, or even larger.
|
|
|
|
The pseudo randomness of numerical random number generators also has
|
|
an advantage. They allow to repeat exactly the same sequence of
|
|
random numbers. After defining the state of the generator or setting a
|
|
\enterm{seed} with the \code{rng()} function, a particular sequence of
|
|
random numbers is generated by subsequent calls of the random number
|
|
generator. This way we can not only precisly define the statistics of
|
|
noise in our simulated data, but we can repeat an experiment with
|
|
exactly the same sequence of noise values. This is useful for plots
|
|
that involve some random numbers but should look the same whenever the
|
|
script is run.
|
|
|
|
\begin{exercise}{randomnumbers.m}{randomnumbers.out}
|
|
First, read the documentation of the \varcode{rand()} function and
|
|
check its output for some (small) input arguments.
|
|
|
|
Generate three times the same sequence of 10 uniformly distributed
|
|
numbers using the \code{rand()} and \code{rng()} functions.
|
|
|
|
Generate 10\,000 uniformly distributed random numbers and compute
|
|
the correlation coefficient between each number and the next one in
|
|
the sequence. This is the serial correlation at lag one.
|
|
\end{exercise}
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Univariate data}
|
|
A basic type of stochastic simulations is to draw random numbers from
|
|
a given distribution like, for example, the normal distribution. This
|
|
simulates repeated measurements of some quantity (weight of tigers,
|
|
firing rate of neurons). Doing so we must specify from which
|
|
probability distribution the data should originate from and what are
|
|
the parameters of that distribution (mean, standard deviation, shape
|
|
parameters). How to illustrate and quantify univariate data, no matter
|
|
whether they have been actually measured or whether they have been
|
|
simulated as described in the following, is described in
|
|
chapter~\ref{descriptivestatisticschapter}.
|
|
|
|
\subsection{Normally distributed data}
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics{normaldata}
|
|
\titlecaption{\label{normaldatafig} Generating normally distributed
|
|
data.}{With the help of a computer the weight of 300 tigers can be
|
|
measured in no time using the \code{randn()} function (left). By
|
|
construction we then even know the population distribution (red
|
|
line, right), its mean (here 220\,kg) and standard deviation
|
|
(40\,kg) from which the simulated data values were drawn (yellow
|
|
histogram).}
|
|
\end{figure}
|
|
|
|
For drawing numbers from a normal distribution
|
|
\begin{equation}
|
|
p_{norm}(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^2}
|
|
\end{equation}
|
|
with zero mean and unit standard deviation we use the \code{randn()}
|
|
function. This function returns normally distributed numbers $\xi_i$
|
|
with zero mean and unit standard deviation. For changing the standard
|
|
deviation we need to multiply the returned data values with the
|
|
required standard deviation $\sigma$. For changing the mean we just
|
|
add the desired mean $\mu$ to the random numbers
|
|
(\figref{normaldatafig}):
|
|
\begin{equation}
|
|
x_i = \sigma \xi_i + \mu
|
|
\end{equation}
|
|
|
|
\begin{exercise}{normaldata.m}{normaldata.out}
|
|
Write a little script that generates $n=100$ normally distributed
|
|
data simulating the weight of Bengal tiger males with mean 220\,kg
|
|
and standard deviation 40\,kg. Check the actual mean and standard
|
|
deviation of the generated data. Do this, let's say, five times
|
|
using a for-loop. Then increase $n$ to 10\,000 and run the code
|
|
again. It is so simple to measure the weight of 10\,000 tigers for
|
|
getting a really good estimate of their mean weight, isn't it?
|
|
Finally plot from the last generated data set of tiger weights the
|
|
first 1000 values as a function of their index.
|
|
\end{exercise}
|
|
|
|
\subsection{Other probability densities}
|
|
We can draw random numbers not only from normal
|
|
distributions. Functions are provided that let you draw random numbers
|
|
of almost any probability distribution. They differ in their shape and
|
|
in the number of parameters they have. Most probability distributions
|
|
are parameterized by a location parameter that usually coincides with
|
|
their mean, and a scale parameter that often coincides with their
|
|
standard deviation. Some distributions have additional shape
|
|
parameters.
|
|
|
|
For example, the time intervals $t$ between randomly generated action
|
|
potentials (a so called Poisson spike train, see
|
|
chapter~\ref{pointprocesseschapter}) are exponentially distributed ---
|
|
as are the intervals between state transitions of ion channels, or the
|
|
intervals between radioactive decays. The exponential distribution is
|
|
only defined for positive time intervals:
|
|
\begin{equation}
|
|
\label{expdistribution}
|
|
p_{exp}(t) = \lambda e^{-\lambda t} \; , \quad t \ge 0, \; \lambda > 0
|
|
\end{equation}
|
|
A single rate parameter $\lambda$ measured in Hertz parameterizes the
|
|
exponential distribution. It defines how often per time unit the event
|
|
happens. Both the mean interval between the events and the
|
|
corresponding standard deviation equal the inverse rate.
|
|
|
|
\begin{exercise}{exprandom.m}{exprandom.out}
|
|
Draw $n=200$ random time intervals from an exponential distribution
|
|
with rate $\lambda=50$\,Hz. Calculate the mean and the standard
|
|
deviation of the random numbers and compare them with the expected
|
|
values.
|
|
\end{exercise}
|
|
|
|
|
|
\subsection{Uniformly distributed random numbers}
|
|
|
|
\begin{figure}[t]
|
|
Exponential distribution of ISIs.
|
|
\end{figure}
|
|
|
|
The \code{rand()} function returns uniformly distributed random
|
|
numbers between zero and one
|
|
\begin{equation}
|
|
p_{uni}(x) = \left\{\begin{array}{rcl} 1 & ; & 0 \le x < 1 \\
|
|
0 & ; & \text{else} \end{array}\right.
|
|
\end{equation}
|
|
For random numbers $x_i$ uniformly distributed between $x_{\rm min}$
|
|
and $x_{\rm max}$ we multiply by the desired range $x_{\rm max} -
|
|
x_{\rm min}$ and add $x_{\rm min}$ to the numbers $u_i$ returned by
|
|
\code{rand()}:
|
|
\begin{equation}
|
|
x_i = (x_{\rm max} - x_{\rm min}) u_i + x_{\rm min}
|
|
\end{equation}
|
|
|
|
\begin{exercise}{unirandom.m}{unirandom.out}
|
|
Simulate the times of $n=10\,000$ action potentials occuring
|
|
randomly within a time interval of 100\,s using the \code{rand()}
|
|
function. Sort the times using \code{sort()}, compute the time
|
|
intervals between the action potentials (interspike intervals, ISIs)
|
|
using the \code{diff()} function, plot a histogram of the interspike
|
|
intervals in milliseconds, and compute the mean and standard
|
|
deviation of the interspike intervals, also in milliseconds. Compare
|
|
the mean ISI with the firing rate (number of action potentials per
|
|
time).
|
|
\end{exercise}
|
|
|
|
|
|
\subsection{Simulating probabilities}
|
|
Simulating events that happen with some probability $P$ is also
|
|
possible. That could be the probability of head showing up when
|
|
flipping a coin, getting an action potential within some time, a ion
|
|
channel to open, ... For this draw a random number $u_i$ from a
|
|
uniform distribution between zero and one (\code{rand()}). If the
|
|
random number is lower than $P$, the event happens, if it is larger
|
|
the event does not occur.
|
|
|
|
\begin{exercise}{probability.m}{probability.out}
|
|
Let's flip a coin 20 times. The coin is biased and shows head with a
|
|
probability of $P=0.6$. Count the number of heads.
|
|
\end{exercise}
|
|
|
|
|
|
The gamma distribution (\code{gamrnd()}) phenomenologically describes
|
|
various types of interspike interval dsitributions
|
|
(chapter~\ref{pointprocesseschapter}). scale and shape. exercise.
|
|
|
|
\subsection{Random integers}
|
|
\code{randi()}
|
|
|
|
|
|
\section{Random walks}
|
|
A basic concept for stochastic models is the \enterm{random walk}. A
|
|
walker starts at some initial position $x_0$. At every iteration $n$
|
|
it takes a step $\delta_n$: $x_n = x_{n-1} + \delta_n$. The walker
|
|
either steps to the right ($\delta_n = +1$) with some probability $p$
|
|
or to the left ($\delta_n = -1$) with probability $q = 1-p$.
|
|
|
|
\begin{figure}[tb]
|
|
\includegraphics{randomwalkone}
|
|
\titlecaption{\label{randomwalkonefig} Random walk in one
|
|
dimension.}{Twelve random walks starting at $x_0=0$. The ensemble
|
|
average of the position of the random walkers marked by the blue
|
|
line and the corresponding standard deviation by the shaded
|
|
area. The standard deviation increases with the square root of the
|
|
number of steps taken. The top plot shows a symmetric random walk
|
|
where the probability of a step to the right equals the one for a
|
|
step to the left ($p=q=\frac{1}{2}$). The bottom plot shows a
|
|
random walk with drift. Here the probability of taking a step to
|
|
the right ($p=0.6$) is larger than the probability for a step to
|
|
the left ($q=0.4$).}
|
|
\end{figure}
|
|
|
|
For a \enterm[random walk!symmetric]{symmetric random walk}
|
|
(\figref{randomwalkonefig} top), the probabilities for a step to the
|
|
left or to the right are the same, $p = q = \frac{1}{2}$. The average
|
|
position of many walkers is then independent of the iteration step and
|
|
equals the initial position $x_0$. The variance of each step is
|
|
$\sigma_{\delta_n}^2 = p \cdot 1^2 + q \cdot(-1)^2 = 1$. Since each
|
|
step is independent of the steps before, the variances of all the
|
|
steps add up such that the total variance of the walker position at
|
|
step $n$ equals $\sigma_{x_n}^2 = \sum_{i=1}^n \sigma_{\delta_n}^2 =
|
|
n$. The standard deviation of the walker positions grows with the
|
|
square root of the number of iteration steps $n$: $\sigma_{x_n} =
|
|
\sqrt{n}$.
|
|
|
|
\enterm[random walk!with drift]{Random walks with drift} are biased to
|
|
one direction (\figref{randomwalkonefig} bottom). The probability to
|
|
step to the left does not equal the probability to step to the right,
|
|
$p \ne q$. Consequently, the average step is no longer zero. Instead
|
|
it equals $\mu_{\delta_n} = p \cdot 1 + q \cdot (-1) = 2p - 1$. The
|
|
average position grows proportional to the number of iterations:
|
|
$\mu_{x_n} = n(2p - 1)$. The variance of each step turns out to be
|
|
$\sigma_{\delta_n}^2 = 4pq$ (this indeed equals one for
|
|
$p=q=\frac{1}{2}$). Again, because of the independence of the steps,
|
|
the variance of the position at the $n$-th iteration, $\sigma_{x_n}^2
|
|
= 4pqn$, is proportional to $n$ and the standard deviation,
|
|
$\sigma_{x_n} = 2\sqrt{pqn}$, grows with the square root of $n$.
|
|
|
|
\begin{exercise}{}{}
|
|
Random walk exercise.
|
|
\end{exercise}
|
|
|
|
\begin{figure}[tb]
|
|
\includegraphics{randomwalkneuron}
|
|
\titlecaption{\label{randomwalkneuronfig} Random walk spiking neuron.}{.}
|
|
\end{figure}
|
|
|
|
The sub-threshold membrane potential of a neuron can be modeled as a
|
|
one-dimensional random walk. The thousands of excitatory and
|
|
inhibitory synaptic inputs of a principal neuron in the brain make the
|
|
membrane potential to randomly step upwards or downwards. Once the
|
|
membrane potential hits a threshold an action potential is generated
|
|
and the membrane potential is reset to the resting potential. The
|
|
threshold in this context is often called an ``absorbing
|
|
boundary''. The time it takes from the initial position to cross the
|
|
threshold is the ``first-passage time''. For neurons this is the
|
|
interspike interval. For symmetric random walks the first-passage
|
|
times are again exponentially distributed. ??? Is that so ???
|
|
|
|
% Gerstein and Mandelbrot 1964
|
|
|
|
\begin{exercise}{}{}
|
|
Random walk neuron exercise.
|
|
\end{exercise}
|
|
|
|
\begin{figure}[tb]
|
|
\includegraphics{randomwalktwo}
|
|
\titlecaption{\label{randomwalktwofig} Random walk in two
|
|
dimensions.}{Fivehundred steps of five random walkers starting at
|
|
the origin (black dot) and taking steps of distance one either up,
|
|
down, to the left, or to the right with equal probability.}
|
|
\end{figure}
|
|
|
|
Higher dimensions (\figref{randomwalktwofig}), Brownian motion. Crystal growth.
|
|
|
|
\begin{exercise}{}{}
|
|
Random walk in 2D exercise.
|
|
\end{exercise}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Bivariate data}
|
|
|
|
\subsection{Static nonlinearities}
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics{staticnonlinearity}
|
|
\titlecaption{\label{staticnonlinearityfig} Generating data
|
|
fluctuating around a function.}{The conductance of the
|
|
mechontransducer channels in hair cells of the inner ear is a
|
|
saturating function of the deflection of their hairs (left, red
|
|
line). Measured data will fluctuate around this function (blue
|
|
dots). Ideally the residuals (yellow histogram) are normally
|
|
distributed (right, red line).}
|
|
\end{figure}
|
|
|
|
Example: mechanotransduciton! (\figref{staticnonlinearityfig})
|
|
|
|
draw (and plot) random functions.
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Dynamical systems}
|
|
|
|
\begin{itemize}
|
|
\item iterated maps
|
|
\item euler forward, odeint
|
|
\item introduce derivatives which are also needed for fitting (move box from there here)
|
|
\item Passive membrane
|
|
\item Add passive membrane to mechanotransduction!
|
|
\item Integrate and fire
|
|
\item Fitzugh-Nagumo
|
|
\item Two coupled neurons? Predator-prey?
|
|
\end{itemize}
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Summary}
|
|
with outook to other simulations (cellular automata, monte carlo, etc.)
|
|
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\printsolutions
|