Methods in Computational Neuroscience

Obidos, 2004

 

 

MATLAB tutorial

 

Assembled, edited and written by Oren Shriki, Oren Farber and Dmitri Bibitchkov.

 

 

 

 

 

 

 

 

 

 

 

 

MATLAB is a general-purpose mathematical software that is user friendly and very useful for data analysis, simulations, presentation of results, and more. The software runs on both UNIX and WINDOWS.

MATLAB is an interpreter, which means that it performs every command line immediately after it is entered. So you can either enter your commands line by line (prompt mode), or prepare a script file first, and then run a program. A script file for MATLAB is called an m-file and it must have the extension 'm' (i.e. something like myprogram.m). To run an m-file you just enter the name of the file without the extension at the MATLAB prompt. Graphical output is available to supplement numerical results.

 

For those who prefer to find out everything themselves!

Online help is available from the Matlab prompt (a double arrow), both generally (listing all available commands):

>> help

[a long list of help topics follows]

 

and for specific commands:

>> help fft

[a help message on the fft function follows in the prompt window]

 

You can also get help on a function or a topic in a separate window using a command helpwin For example,

>> helpwin datafun

opens a desktop help window with information about functions for data analysis and Fourier transforms.  At the end of each text concerning a function you can find links to related functions, indicated by 'See Also'. For general help, for example if you would like to get a general introduction into a specific Toolbox, use the command

>> helpdesk

 

This will display the start page for the online documentation in a MATLAB help browser or in your web browser, depending on the settings of your system. The online documentation is also available from http://www.mathworks.com/access/helpdesk/help/helpdesk.html

 

If you want to search for commands related to some topic, say, complex numbers, you enter

>> lookfor complex

 

 

 

Getting started

Defining matrices

 

 

If you are a first-time user, you should invoke MATLAB now and follow along. In MATLAB you don't have to define a variable type - you just enter it. The main object of the MATLAB is the matrix (The name MATLAB stands for "matrix laboratory").

You can enter matrices into MATLAB in several different ways.

·       Enter an explicit list of elements.

·       Load matrices from external data files.

·       Generate matrices using built-in functions.

·       Create matrices with your own functions in M-files.

 

A matrix of one element is called a scalar and we can create one by entering:

>> s = 3

 

If you enter a variable this way, MATLAB immediately displays its value in the prompt window. To avoid this, put the ’;’ at the end of the line. You can use the ';' also to separate different commands on the same line like this

>> s = 3; a=10;

 

A matrix of one row or one column is called a vector. We can create a row vector by entering:

>> v = [8 3 5[

The result is a 1 by 3 matrix

 

To create a column vector we will use the ';' symbol to separate between the elements

>> u = [4.2; 5; 7]

The result is a 3 by 1 matrix

 

We can create a 3 by 3 matrix simply by entering

>> A = [3 4 0; 1 4 6; 9 7 2]

 

We can now calculate quantities like the transpose of A

>> A'

or the transpose of v which will be a column vector:

>> v'

 

For example we can create another 3 by 3 matrix by entering

>> B=[v; u'; 1 0 4[

 

Note that if A has complex entries then A' means a complex conjugate transpose. A non-conjugate transpose is performed by

>> A.'

Check this by typing

[X, Y]=meshgrid(-1:1);

Z=X+i*Y

Z.'-Z

 

MATLAB has some built-in functions for the construction of special matrices.

 

>> ones(3,4)

will create a matrix of 1's with 3 rows and 4 columns. And you can figure out what

>> zeros(5,7(

does.  

 

As you see, every Matlab function (even those you will write) has input arguments and output arguments . The input arguments are what you put between the brackets, for example in ones(m,n) m and n are the inputs. The output is what you get back from the function. For example in A = ones(5) the result of the ones function is placed in A , which will be  a 5-by-5 matrix of zeros. The variables you use to get the output are the glue between the various commands. In order to find about the input and output of functions, simply type the magic word help.

 

A command

>>rand(3,4)

 

creates a 3 by 4 matrix with random values chosen from a uniform distribution on the interval (0.0,1.0).

 

>>repmat(a,M,N)

creates an M by N matrix filled by values of a. It produces a large matrix consisting of M by N copies of a, if  a is a matrix. To get an array of copies of a matrix A, try

>>B=repmat(A,[1 1 4])

 

and see what it does.

 

To create an identity matrix of order 5 you type

>> eye(5)

 

Another  example matrix appears in the Renaissance engraving Melancholia I by the German artist and amateur mathematician Albrecht Dürer. This matrix is known as a magic square and has a property  that it has equal row, column, and diagonal sums. You can create such a matrix using a function magic

>> A=magic(4)

 

creates an 4-by-4 magic square.

 

If you type

>> x=1:10

this will give a vector with all the integer numbers from 1 to 10, and

>> y=0:0.1:8

will give a vector of all the numbers between 0 and 8 with a 0.1 step.

 

To summarize all we have learned, enter

>> s=3; y=[1 s ones(1,5); 5*rand(1,7); 1:7[;

 

and type y to see what you've got.

 

Indexing

 

After we created a matrix we might want to extract some elements from it. To extract the j'th element of a vector you type

>> v(j)

Typing the command  A(j,k)  will extract the element of A that sits in the j'th row in the k'th column. For example:

>> A(2,3)

will extract the element in the 2nd row in the 3rd column,  and

>> A(:,3:end(

will extract the last columns of A beginning from the 3rd.  The first semicolon stands for all the indexes along the corresponding dimension and end denotes the last index in an indexing expression.

 

It is also possible to refer to the elements of a matrix with a single subscript, A(k). This is the usual way of referencing row and column vectors. But it can also apply to a fully two-dimensional matrix, in which case the array is regarded as one long column vector formed from the columns of the original matrix. So, for our magic square, A(8) is another way of referring to the value 15 stored in A(4,2). If you try to use the value of an element outside of the matrix, it is an error:

>> t = A(4,5)

Index exceeds matrix dimensions.

On the other hand, if you store a value in an element outside of the matrix, the size increases to accommodate the newcomer:

>> X = A;

>> X(4,5) = 17

 

X =

    16     3     2    13     0

     5    10    11     8     0

     9     6     7    12     0

     4    15    14     1    17

 

Operations with matrices

 

If we have several matrices, we can add them

>> A+B

or multiply them

>> A*B

 

but we have to take care that the dimensions are appropriate: the number of columns in A should match the number of rows in B. The multiplication operation we just used obeys the rules of ordinary matrix multiplication. There is another kind of multiplication, which is multiplying each element in a matrix by the corresponding element of  another matrix which has to be the same size (this is called array multiplication). To do this you enter

>> A.*B

 

This dot rule is general, so if you enter

>> A.^2

 

it will raise all the elements of A to the power of 2 (which is of course different from A^2).

MATLAB uses a dot, or decimal point, as part of the notation for element-by-element operations.

The list of operators includes:

+   

Addition

-   

Subtraction

.*

Element-by-element multiplication

./

Element-by-element division

.\

Element-by-element left division

.^

Element-by-element power

.'  

Unconjugated array transpose

 

A command

>>sum(A)

performs a summation along columns. Suppose we want to know the sum of the 3rd row. We will use the colon (:) and the function sum like this:

>> sum(A(3,:))

for the sum of the third column we would type:

>> sum(A(:,3))

 

The main diagonal of the matrix can be picked off by a function diag and the determinant is calculated by function det. The other diagonal, the so-called antidiagonal, is not so important mathematically, so MATLAB does not have a ready-made function for it. But a function originally intended for use in graphics, fliplr flips a matrix from left to right.

>>sum(diag(fliplr(A)))

 

calculates the sum of the antidiagonal

 

Now, type the following command:

>> A>3

What you got is a new matrix, with ones and zeros. Every element with one means that the corresponding element in A fulfills the condition "bigger than 3". Zeros mean that the opposite is true, that is, the number is smaller or equal to 3. We now go further by typing:

>> (A>3).* A

can you explain the result ? 

hint: notice we used the element by element multiplier (.*)

 

Now try this and explain the result:

>> A.*A > 5*A

 

Logical Subscripting

 

The logical vectors created from logical and relational operations can be used to reference subarrays. Suppose X is an ordinary matrix and L is a matrix of the same size that is the result of some logical operation. Then X(L) specifies the elements of X where the elements of L are nonzero.

This kind of subscripting can be done in one step by specifying the logical operation as the subscripting expression. Suppose you have the following set of data.

x =

  2.1 1.7 1.6 1.5 NaN 1.9 1.8 1.5 5.1 1.8 1.4 2.2 1.6 1.8

 

The NaN is a marker for a missing observation, such as a failure to respond to an item on

a questionnaire. To remove the missing data with logical indexing, use finite(x), which is true for all finite numerical values and false for NaN and Inf.

x = x(finite(x))

x =

  2.1 1.7 1.6 1.5 1.9 1.8 1.5 5.1 1.8 1.4 2.2 1.6 1.8

 

Now there is one observation, 5.1, which seems to be very different from the others. It is an outlier. The following statement removes outliers, in this case those elements more than three standard deviations from the mean.

x = x(abs(x-mean(x)) <= 3*std(x))

x =

  2.1 1.7 1.6 1.5 1.9 1.8 1.5 1.8 1.4 2.2 1.6 1.8

 

The find Function

 

The find function determines the indices of array elements that meet a given logical condition. In its simplest form, find returns a column vector of indices. Transpose that vector to obtain a row vector of indices. For example

k = find(isprime(A))'

 

picks out the locations, using one-dimensional indexing, of the primes in the magic square.

Display those primes, as a row vector in the order determined by k, with

A(k)

 

When you use k as a left-hand-side index in an assignment statement, the matrix structure is preserved. For example,

A(k) = NaN

 

will change previously prime entries of A into NaN, with the rest of the matrix left unchanged.

Data storage

 

While working with MATLAB, you can always list the variables in the workspace by entering who or whos ( whos shows also the size of the variables)

To log everything you type in the workspace (and the results), use the 'diary' command. To start, enter something like:

>> diary my_log.txt

 

and to stop saving, enter:

>> diary off

 

Type help diary to see exactly what it does.

 

You can store all your data currently defined in the workspace in a file :

>> save mydata

 

a file mydata.mat will be created with  all the variables you currently have. The extension .mat   indicates the MATLAB data format. Note, that if this file already exists, MATLAB will overwrite the old data. If you want to add some variable to an existing file mydata.mat type

>>save mydata A -append

 

You can remove all the defined variables using the command clear and then load all the data again from file :

>> clear

>> load mydata

 

You can also load and store data using the MATLAB user interface. To save all variables in the workspace, choose a corresponding option in the File menu. If you want to save a certain variable, choose it in the Workspace window, click the right mouse button and select the appropriate option.

You can load a data file by selecting it in the Current Directory window. However, this does not allow you to load just selected variables from the file, which is possible using the command load:

>> load mydata A x

 

It is also possible to load the data stored in other formats rather than MATLAB. For this you need to Import data using the Import Wizard dialog. You can start a demo Importing Data from the helpdesk to get more information about this tool.

 

Data Analysis

MATLAB uses column-oriented analysis for multivariate statistical data. It assumes that each column in a data set represents a variable and each row an observation. The (i,j)th element is the ith observation of the jth variable. For vector arguments, it does not matter whether the vectors are oriented in a row or column direction. For array arguments, the functions operate on columns. This means, for example, that if you apply max to an array, the result is a row vector containing the maximum values over each column.

 

 

Here we list few functions, which are very useful for data analysis:

max - maximum value

min - minimum value

mean - mean value

median - median value

std - standard deviation

sort – sorting

sum -sum of the elements

prod - product of elements

cumsum - cumulative sum of elements

cumprod - cumulative product of elements

diff - differences and approximate derivatives

corrcoef - correlation coefficients

cov - covariance matrix

 

Typing help on each of these functions will give you the full description of them. Note that, although it is very easy to write most of these functions yourself, it is nice to already have them and it makes the scripts more compact.

As an example, consider a data set with three variables:

·       Heart rate

·       Weight

·       Hours of exercise per week

For five observations, the resulting array might look like:

D =

       72          134          3.2

       81          201          3.5

       69          156          7.1

       82          148          2.4

       75          170          1.2

To obtain the mean and standard deviation of each column:

>> mu = mean(D), sigma = std(D)

 

For a list of the data analysis functions available in MATLAB, type

>> help datafun

If you have access to the Statistics Toolbox, type

>> help stats

 

To learn about the interpolation capabilities of MATLAB, try the following demo

>> census

 

 

 

Simple statistics

Hypothesis testing for the difference in means of two samples.

Syntax

[h,significance,ci] = ttest2(x,y)

[h,significance,ci] = ttest2(x,y,alpha)

[h,significance,ci] = ttest2(x,y,alpha,tail)

Description

h = ttest2(x,y) performs a t-test to determine whether two samples from a normal distribution (in x and y) could have the same mean when the standard deviations are unknown but assumed equal.

h, the result, is 1 if you can reject the null hypothesis at the 0.05 significance level alpha and 0 otherwise.

significance is the p-value associated with the T-statistic.

significance is the probability that the observed value of T could be as large or larger under the null hypothesis that the mean of x is equal to the mean of y.

ci is a 95% confidence interval for the true difference in means.

[h,significance,ci] = ttest2(x,y,alpha) gives control of the significance level, alpha. For example if alpha = 0.01, and the result, h, is 1, you can reject the null hypothesis at the significance level 0.01. ci in this case is a
100(1-alpha)% confidence interval for the true difference in means.

ttest2(x,y,alpha,tail) allows specification of one or two-tailed tests. tail is a flag that specifies one of three alternative hypotheses:

·       tail = 0 (default) specifies the alternative, .

·       tail = 1 specifies the alternative, .

·       tail = -1 specifies the alternative, .

Examples

This example generates 100 normal random numbers with theoretical mean zero and standard deviation one. We then generate 100 more normal random numbers with theoretical mean one half and standard deviation one. The observed means and standard deviations are different from their theoretical values, of course. We test the hypothesis that there is no true difference between the two means. Notice that the true difference is only one half of the standard deviation of the individual observations, so we are trying to detect a signal that is only one half the size of the inherent noise in the process.

x = normrnd(0,1,100,1);

y = normrnd(0.5,1,100,1);

[h,significance,ci] = ttest2(x,y)

 

h =

1

 

significance =

0.0017

 

ci =

-0.7352 -0.1720

 

The result, h =1, means that we can reject the null hypothesis. The
significance is 0.0017, which means that by chance we would have
observed values of t more extreme than the one in this example in only 17
of 10,000 similar experiments! A 95% confidence interval on the mean is
[-0.7352 -0.1720], which includes the theoretical (and hypothesized) difference of -0.5.

 

Graphics and Visualization

 

MATLAB is indeed a great tool for visualization of data of all sorts. The first command you should learn how to use is the plot. If you want to plot vector y versus vector x just type

>> plot(x,y(

MATLAB will produce a figure window and will plot the graph.

For example, type the following:

>> x = 0:pi/8:pi*2;

>> m=sin(x);

>> plot(x, m)

 

To learn more about the plot command, try the help on it. To present several plots on the same figure, use the subplot command.  To create another figure window, use the figure command.

Now suppose we have a  normally distributed random variable y with the mean, which depends on the parameter x as sin(x), and has a constant standard deviation 0.1. Let's produce 100 'observations'  of  this variable per each x:

>>y=repmat(sin(x),100,1)+randn(100,length(x))*0.1;

 

Now, we plot the data:

>> plot(x,y,'.');

 

The 3rd argument specifies the symbol with which we plot the data. See help to find out what are the other options. Now, let's create a plot of the data with error bars:

>> errorbar(x,mean(y),std(y));

 

Now, we shall learn about a useful command: hist.  With  hist you get nice histograms very quickly. Try this:

>> x = rand(1000,1) ;

>> hist(x)

 

this will create a histogram of a random vector with 1000 elements.

>> max(x)

this is the greatest element

>> max(x) – min(x)

this is the difference between the greatest and lowest. A command range produces the same result. We can define the binning ourselves the following way:

>> dx=min(x):range(x)/100:max(x);

>> h=hist(x,dx);

>> bar(dx,h);

The command bar creates a bar graph like the one produced by the command hist.

 

The following commands will help you control the way your graphs look: axis, xlabel, ylabel, title, hold, grid, get, set, gca, gcf. You can also edit your graphs manually and save figures in various formats using the interactive menu of the window, in which your figure is displayed. An interactive tool called Properties Editor enables you to see and change object property values.

A 3-D graph is plotted using commands surf and mesh. To illustrate this, consider the following example. The strength of connections between columns of neurons is often assumed to be dependent on a distance between these columns in form of a so-called Mexican-hat function, characterized by local excitation and long-range inhibition.  Let us demonstrate such interactions of a neuron placed at the origin with its surrounding on an area of 3-by-3 mm

 

x=-3:0.05:3;

[X, Y]=meshgrid(x);

R=X.^2+Y.^2;

surf(X,Y,R);

pause

surf(X,Y,exp(-R/0.5)-0.5*exp(-R/3));

 

A command meshgrid returns the X- and Y- coordinates of the nodes of a grid built on the domain specified by vector x.    The rows of the output array X and  the columns of the output array Y  are copies of the vector x. You can rotate the graph by clicking the corresponding button in the figure menu and then dragging the graph with the mouse. This way you can see the graph from any angle, especially you can observe any 2-D projection.

 

The best way to learn about MATLAB visualization capabilities is to see the demo and then learn about the commands that are most useful for you. Note that the demos are written in MATLAB so you can "type" the files and just copy the commands you need.

 

Functions and scripts

Scripts

 

Instead of entering all the commands in the command line you can write them all in an M-file. You create MATLAB scripts using a text editor, then use them as you would any other MATLAB  command. Suppose you have created a script  myscript.m , then you can simply run it from the command  prompt by typing

 

>> myscript

 

or use it in other scripts. For this,  you need to run both scripts from the same directory or to specify the path to the directory where  myscript.m is stored. Suppose you have stored your files in the directory C:\MATLAB6p5\work. To be able to access these files from other directories, type

>>addpath C:\MATLAB6p5\work

or set the path to the directory using the File menu.

 Scripts do not accept input arguments or return output arguments. They operate on data in the workspace or can create new data on which to operate. Any variables that they create remain in the workspace, to be used in subsequent computations.

 

Functions

 

You can define your own functions which can accept input arguments and return output arguments.

 Suppose, you have just discovered the formula that explains the universe and everything:

Now you can create a file called universe.m with the following content

function y=universe(x);

y=x.^2+42*sin(x);

 

Note, that the name of the function must be identical with the name of the file! The variable y defines the output variable of the function. This is what you get back from the function. Unlike in a script, the rest of the variables used by the function are internal and are not passed to the workspace from which the function was called.  The input arguments are defined between the brackets.  Your function can have several input arguments. For example, you can define an optional parameter a and set its default value to 42:

function y=universe1(x,a);

if nargin < 2

        a=42;

end

y=x.^2+a*sin(x);

 

Inside the body of a user-defined function, the variable nargin counts the number of parameters actually used to call the function. Now, if you call the function universe(x) it will not return an error but will calculate the value using a=42. Another way to check if a certain variable is defined is to use a function exist:

if ~exist('a','var')

        a=42;

end

 

Finally, you can declare a function with variable number of arguments using varargin.

 

function y=universe2(x,varargin)

 

varargin must be declared as the last input argument    and collects all the inputs from that point onwards. Then we can use them the following way:

a=42;

b=1;

if nargin ==2

        a=varargin{1};

elseif nargin>2

      a=varargin{1};

        b=varargin{2};

end

y=b*x.^2+a*sin(x)

 

 

Studying function properties

 

Type:  fplot(‘x^2+42*sin(x)’, [-3*pi 3*pi]).  As you guessed, this plots a function whose definition is specified by the first argument (a string). We  could also  enter 'universe' as a first argument, this would plot our function with its default values for a and b. The second argument is the x range of interest.  Type help fplot to learn more.

Try to plot the function over larger and smaller ranges. Try to use the zoom tool to inspect the function in details. Add grids to the plot by : grid on (you can turn it of by grid off).

There is another way to define our function, which is especially suitable if we want to use the function as an argument of other functions of MATLAB: 

  f = inline('x^2 + 42*sin(x)')

 

You will get this answer:       

f =

     Inline function:

     f(x) = x^2+42*sin(x)

 

Which means MATLAB now recognizes our function. In order to know how much the universe will last, you need to know for what x between –2 and –4 the function zeros. Type

t=fzero (f, [-2 –4]))

 

If you want to test this, simply type: f(t). Use this function with different second arguments to find other intersections with zero. Now feel a strong urge to find the minimum points of f. This is simple.

Type: m=fmin(f, -8, 0)

This means that m will be assigned with the x that gives a local minimum such that   -8<m<0.

In order to perform analytical calculations of any kind, one can use a Symbolic Math Toolbox.

For example, we can calculate the derivative of f. Type:

s = diff('x^2 + 42*sin(x)');

 

The output will be a symbolic expression containing the derivative of the input, if the argument of the function is a string, or a symbolic expression (but not an inline function).  You can define a symbolic variable by typing

>> syms x

 

Then we don't need to use '':

 

s = diff(x^2 + 42*sin(x));

 

 If we want to plot the result, we need to use the string equivalent to s:

ss = char(s)

figure    (this will create another figure)

fplot(ss, [-3*pi 3*pi])

 

If we wanted to plot the integral f we would do

s = int(x^2 + 42*sin(x))

 

and repeat the other command as before.

 

Want it faster ?  C, seňor !

Speed of simulations is usually considered as disadvantage of MATLAB. It is especially slow if it has to operate with  for - loops and allocate large amount of memory. One of the possibilities to make programs work faster and still enjoy the advantages of MATLAB interface is to write time- and memory-consuming parts of the code in compiled languages like Fortran and C and use them as functions in MATLAB scripts. MATLAB gives this possibility via MEX-files. MEX-files are dynamically linked subroutines produced from C or Fortran source code that the MATLAB interpreter can automatically load and execute. They behave just like M-files and built-in functions. While M-files have a platform-independent extension, .m, MATLAB identifies MEX-files by platform-specific extensions. For Windows, for example, it is .dll (dynamically linked libraries).

Before you can create MEX-files on the Windows platform, you must define which compiler you are going to use. For this, type

mex -setup

 

You will get a list of compilers installed on your platform.  MATLAB includes a C compiler, Lcc, that you can use to create C MEX-files if you don't have any other compiler installed.

Next, you need to create a C-code. The C code for a MEX-file differs from usual C programs, because it should be able to convert MATLAB variables passed to the function into types "digestable" by C, and transform the result of computations back into a MATLAB array.

The source code for a MEX-file consists of two distinct parts: A computational routine that contains the code for performing the computations that you want implemented in the MEX-file. Computations can be numerical computations as well as inputting and outputting data. A gateway routine that interfaces the computational routine with MATLAB by the entry point mexFunction (see example). Its parameter right is an array of  input arguments, n_right is the number of  input arguments, left is an array of output arguments, and n_left is the number of output arguments. The gateway calls the computational routine as a subroutine.

Here is an example of a function which acts like square root for positive inputs, but returns negative values for negative inputs instead of imaginary. It is actually equivalent to sign(x).*sqrt(abs(x));

 

#include <stdlib.h>

#include <math.h>

#include "mex.h"

void negsqrt(double *y, double *x, int n)

{    

        int i;

        double dummy;

        for (i=0; i<n; ++i)

        {

                 dummy = x[i];

                 y[i] = ((dummy < 0) ? -sqrt(-dummy) : sqrt(dummy));

        }

        return;

}

 

void mexFunction(int n_left, mxArray *left[],

        int n_right, const mxArray *right[])

{

        int  n;

        double *source, *destination, dummy;

 

        if (n_left != 1)

                 mexErrMsgTxt("Please enter the output array!\n");

        if (n_right != 1)

                 mexErrMsgTxt("Please enter the input array!\n");

 

        if (!mxIsDouble(right[0]))

                 mexErrMsgTxt("Input should be double!\n");

 

        n = mxGetNumberOfElements(right[0]);

        left[0] = mxCreateDoubleMatrix(1, n, mxREAL);

        source = mxGetPr(right[0]);

        destination = mxGetPr(left[0]);

       

/*      do something...  */

      negsqrt(destination, source, n);

        return;

}

 

 

                                                                     Exercises

Level 1

Ex1 : Basics

1.     Create a matrix of 100 random pairs in the range [-1, 1]

2.     Plot the numbers with red dots

3.     Calculate the distances vector

4.     Make a histogram of the distances vector

5.     Find the maximal distance

6.     Find the index of the maximal distance

7.     Find the minimal distance

8.     Find the index of the minimal distance

9.     Plot the minimal distance and the maximal distance in blue (don’t forget to type hold on when adding new items to an existing figure)

 

Ex2: Matrix manipulations

1.     Create a matrix with 10 rows and 4 columns. The numbers in the matrix should be random integers within 10 to 100.

2.     Find the maximal values in each column.

3.     Find the maximal values in each row.

4.     Find the maximal entry in the matrix

5.     Find the index of this entry (row and column)

6.     Create a matrix with the same dimensions in which all entries equal the maximal value that you found in 4.

7.     Calculate the difference between the matrix from 6 and the original matrix.

8.     Find all the entries that are less than 40 and more than 20. Can you do it in one line?

 

Ex3: Basic statistics

1.     Create a vector x with 100 random entries between 0 to 1

2.     Create a vector y with 100 random entries between 0.1 to 1.1

3.     Create a figure with two subplots. The upper one should show an histogram of x and the lower an histogram of y. Make sure that both histograms show the same range.

4.     Calculate the mean, range and standard deviation for each vector.

5.     Check the null hypothesis that the two vectors are from the same distribution (use ttest2 or ttest).

6.     Find how many members of x are greater than their corresponding members in y.

7.     Find the greatest difference between x and y.

 

 

Level 2

[P.Dayan and L.F.Abbott, "Theoretical Neuroscience"]

 

1. Generate spikes for 100 s using a Poisson spike generator with a constant rate of 100 Hz, and record their times of occurrence. Compute the coefficient of variation of the interspike intervals, and the Fano factor for spike counts obtained over counting intervals ranging from 1 to 100 ms. Plot the interspike interval histogram.

 

Hints:

For a Poisson process with a firing rate r , the probability of an interspike interval falling between τ and τ +t is

P[τ ≤ ti+1 -  ti < τ +Δt] = rΔt exp(rτ).  Compute spike times ti for i = 1, 2,... iteratively by generating inter-spike intervals from an exponential probability density. If xrand is uniformly distributed over the range between zero and one, the negative of its logarithm is exponentially distributed. Thus, we can generate spike times iteratively from the

formula ti+1 = ti - ln(xrand)/r.

 

The coefficient of variation CV  is the ratio of the standard deviation to the mean. The Fano factor for spike counts  is the ratio between the variance and mean of the spike count:

 

 

2. Add a refractory period to the Poisson spike generator by allowing the firing rate to depend on time. Initially, set the firing rate to a constant value r(t)=r0 . After every spike, set r(t)=0  and then allow it to recover exponentially back r0 with a time constant τref that controls the refractory recovery rate. In other words, have r(t) obey the equation

 

 

except immediately after a spike, when it is set to 0. Plot the coefficient of variation as a function of τref over the range from 1ms to 20ms and plot interspike interval histograms for a few different values of τref. Compute the Fano factor for spike counts obtained over counting intervals ranging from 1 to 100 ms for the case τref =10 ms.

 

Hints: The simplest way to simulate a differential equation is according to Euler approximation. Divide the time into small (compared to 1/r and to τref) time intervals and calculate the rate iteratively:

However, for this simple equation, it is possible to calculate the rate exactly, by integration from the time of the last spike tlast :

 

The algorithm discussed in the exercise 1 only works for constant firing rates. There are two ways to generate trains with a variable firing rate. A simple procedure for generating spikes, which may be more suitable for this exercise, is based on the fact that the estimated probability of firing a spike during a short interval of duration Δt is r(t)·t. The program progresses through time in small steps of size Δt and generates, at each time step, a random number xrand chosen uniformly in the range between zero and one. If r(t)· Δt > xrand at that time step, a spike is fired, otherwise it is not. For the second method, see exercise 3.  

 

 

 

3. Compute autocorrelation histograms of spike trains generated by  Poisson generator with a constant firing rate of 100 Hz, a constant firing rate of 100 Hz together with a refractory period modeled as in exercise 2 with τref=10ms, and a variable firing rate r(t)=100(1+cos(2πt/50ms)).  Plot the histograms over a range from 0 to 100 ms. Approximate the firing rates from these spike trains  in the following ways:

a) by binning time and counting spikes with Δt = 100 ms.

b) using a sliding a rectangular window function along the spike train with Δt = 100 ms.

c) using a Gaussian window function with σt = 100 ms.

d) Approximate firing rate for an α function  window w(τ) = [α2τ exp(-ατ)]+ with 1/α = 100 ms. Hints:

The method used in exercise 1 can be extended to time-dependent rates by using a procedure called the rejection sampling or spike thinning. The thinning technique requires a bound rmax on the estimated firing rate such that r ≤ rmax all times. We first generate a spike sequence corresponding to the constant rate rmax by iterating the rule ti+1 = ti - ln(xrand)/r. The spikes are then thinned by generating another yrand for       each i and removing the spike at time ti from the train if  r(ti)/rmax < yrand. If r(ti)/rmax ≥ yrand, spike i is retained. Thinning corrects for the difference between the estimated time-dependent rate and the maximum rate.

 

 

The spike-train autocorrelation function is defined as

, with

 

The autocorrelation function is just a histogram of all the spike intervals that are present in the data. It is a particular case of the crosscorelation function. Here is a simple MATLAB function for crosscorrelation:

 

function crosscorr(t1,t2,N_bins);

 

M1 = length(t1);

M2 = length(t2);

D = ones(M2,1)*t1 - t2'*ones(1,M1);

D = D(:);

hist(D,N_bins)

 

For rate approximation, use function conv which calculates a convolution of your spike train with a vector containing the values of the window function on the appropriate time interval. The central element of the vector should correspond the time delay τ=0.