diff --git a/simulations/code/exprandom.m b/simulations/code/exprandom.m index 6e20a9b..be9c5a3 100644 --- a/simulations/code/exprandom.m +++ b/simulations/code/exprandom.m @@ -1,8 +1,8 @@ -n = 10000; -lambda = 50.0; % event rate -t = exprnd(1.0/lambda, n, 1); % exponentially distributed random numbers -m = mean(t); % mean time interval -s = std(t); % corresponding standard deviation -fprintf('Generated n=%d exponentially distributed random time intervals\n', n); -fprintf(' with rate=%.0fHz (-> mean=%.3fs, std=%.3fs)\n', lambda, 1.0/lambda, 1.0/lambda); -fprintf('Measured mean=%.3fs, std=%.3fs\n', m, s); +n = 200; % number of time intervals +rate = 50.0; % event rate +t = exprnd(1.0/rate, n, 1); % exponentially distributed intervals +m = mean(t); % mean time interval +s = std(t); % corresponding standard deviation +fprintf('Generated n=%d random time intervals with rate=%.0fHz:\n', n, rate); +fprintf(' Expected mean interval=%.3fs, std=%.3fs\n', 1.0/rate, 1.0/rate); +fprintf(' Measured mean interval=%.3fs, std=%.3fs\n', m, s); diff --git a/simulations/code/exprandom.out b/simulations/code/exprandom.out index 2c91625..7053ee8 100644 --- a/simulations/code/exprandom.out +++ b/simulations/code/exprandom.out @@ -1,3 +1,3 @@ -Generated n=10000 exponentially distributed random time intervals - with rate=50Hz (-> mean=0.020s, std=0.020s) -Measured mean=0.020s, std=0.020s +Generated n=200 random time intervals with rate=50Hz: + Expected mean interval=0.020s, std=0.020s + Measured mean interval=0.020s, std=0.019s diff --git a/simulations/code/unirandom.m b/simulations/code/unirandom.m new file mode 100644 index 0000000..0c6f1fb --- /dev/null +++ b/simulations/code/unirandom.m @@ -0,0 +1,16 @@ +n = 10000; % number of action potentials +T = 100.0; % total time interval +spikes = sort(T*rand(n, 1)); % sorted times of action potentials +isis = diff(spikes); % interspike intervals +rate = n/T; % firing rate +misi = mean(isis); % mean interspike interval +sisi = std(isis); % and standard deviation +fprintf('firing rate = %.1fHz\n', rate); +fprintf(' mean ISI = %.1fms\n', 1000.0*misi); % inverse of rate +fprintf(' std ISI = %.1fms\n', 1000.0*sisi); % same as mean +hist(1000.0*isis, 50); % exponential distribution +xlabel('ISI [ms]'); +ylabel('count'); + + + diff --git a/simulations/code/unirandom.out b/simulations/code/unirandom.out new file mode 100644 index 0000000..4abe24c --- /dev/null +++ b/simulations/code/unirandom.out @@ -0,0 +1,3 @@ +firing rate = 100.0Hz + mean ISI = 10.0ms + std ISI = 10.0ms diff --git a/simulations/lecture/simulations.tex b/simulations/lecture/simulations.tex index 15cb87f..c4dbfa1 100644 --- a/simulations/lecture/simulations.tex +++ b/simulations/lecture/simulations.tex @@ -79,24 +79,28 @@ script is run. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Univariate data} -The most basic type of simulation is to draw random numbers from a -given distribution like, for example, the normal distribution. This -simulates repeated measurements of some quantity (e.g., weight of -tigers or firing rate of neurons). Doing so we must specify from which +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 +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} -For drawing numbers $x_i$ from a normal distribution 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: +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: \begin{equation} x_i = \sigma \xi_i + \mu \end{equation} @@ -144,24 +148,44 @@ only defined for positive time intervals: \label{expdistribution} p_{exp}(t) = \lambda e^{-\lambda t} \; , \quad t \ge 0, \; \lambda > 0 \end{equation} -The exponential distribution is parameterized by a single rate -parameter $\lambda$ measured in Hertz. 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. +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=10\,000$ 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. + 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} -The gamma distribution (\code{gamrnd()}) phenomenologically describes -various types of interspike interval dsitributions -(chapter~\ref{pointprocesseschapter}). scale and shape. exercise. +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} - \code{rand()} between xmin -and xmax. \subsection{Simulating probabilities} Simulating events that happen with some probability $P$ is also @@ -177,6 +201,18 @@ the event does not occur. probability of $P=0.6$. Count the number of heads. \end{exercise} +Random walk! Have a figure with a few 1D random walks and the +increasing standard deviation as a sheded area behind. + + + + + + +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()}