MIT Department of Brain and Cognitive Sciences 9.29J, Spring 2004 -Introduction to Computational Neuroscience Instructor: Professor Sebastian Seung MATLAB, Statistics, and Linear Regression Justin Werfel 9.29 Optional Lecture #1, 2/09/04 1 MATLAB MATLAB is a powerful software package for matrix manipulation. It’s a very useful language not only for this class, but for a variety of scientific applications, and is used widely thoughout industry. Just as when you have a hammer, everything looks like a nail, so when you have MATLAB, everything looks like a matrix. You may learn to approach this Zen-like state by the end of the class. 1.1 Free your mind The basic MATLAB data type is a matrix, an array of values (by default, double-precision floating point numbers), typically arranged as a two-dimensional1 grid (though arrangements of more or fewer dimensions are possible). The simplest matrix consists of a single element. We can assign a value to a variable with a statement like2 i = 4Notice that when you enter this command, MATLAB echoes the result back to you. You can suppress that echo by ending the statement with a semicolon: i = 4;A vector (one-dimensional matrix) can be created with the colon operator. To make a row vector of successive integers, use a statement like x=2:9. The default increment between successive entries is 1; you can specify a different increment by putting it betwwee two colons, e.g., x=9:-2:3 creates a row vector with the elements (9, 7, 5, 3). The transpose operator in MATLAB is the single-quote; the easiest way to create a column vector, then, is to transpose a row vector, e.g., x=(1:8)’. Note that parentheese in MATLAB, like in C and other languages, can be used to specify the order in 1There are at least two ways we can use the word ‘dimension’ when talking about matrices. One refers to the layout, in the sense that a line is one-dimensional, a square two-dimensional, etc. The other is in the sense of dimensionality of a matrix, e.g., (x1 ,x2 ,x3,x4) is four-dimensional (although as a vector, its elements are arranged in a line when it’s written; that is, a vector is one-dimensional in the first sense). MATLAB’s help system tends to use the first sense, Sebastian tends to use the second. Context should hopefully always make it clear which sense is meant. If you have questions about any specific situation, ask. 2I’m not going to display the results of every operation here; try them out for yourself and see what happens. 1 which you want operations performed; if you type x=1:8’, MATLAB first transposes the element 3 (with no effect), then applies the colon operator and gives you a row vector as before. To enter a two-dimensional matrix, use commas or spaces to separate entries on the same row, and semicolons or carriage returns to separate rows. For example, the statement A = [1, 5,3 62 7 4,4; 3, 8 3 8]specifies the matrix ⎣� 1536 74�24� 3838 Parentheses are used to index matrices, and commas to separate dimensions; so for instance, with the last definition we used for x, x(5) is 5, and A(1,4) is 6.3 You can specify ranges of indices with the colon operator; a colon by itself means all entries in that dimension; the keyword end specifies the last entry in a dimension. See what these statements give you: A(2:3,1:3)A(:,2)A(3,2:end)A(1:end-1,3)Some operators have special meanings when applied to matrices (notably multiplicaation *, and raising to a power, ˆ); if you want to apply an operation to eveer element of a matrix separately, put a period before the operator. For instance, (1:8)*(1:8) gives an error because the inner dimensions of the matrices don’t match; (1:8)*(1:8)’is an inner product, returning a single value; (1:8).*(1:8) gives you a row vector of the squares of the first eight counting numbers; so does (1:8).ˆ2. 1.2 M-files, data sets, looping and why you should never do it Suppose we’ve got some set of N sensors from which we can read measurements (temperaature voltage, whatever) all at the same time; and we sample from them repeatedly M times, so we have M sets of N values. Here’s our first look at MATLAB as hammeer one natural way to store these data points is in a matrix A, say M rows by N columns. Then each row is a snapshot of the sensor readings at a given time; A(i, j) is the value the jth sensor had at time i. Let’s further suppose that we want to calculate some weighted sum of the sensor readings at each time step; the weights are given by values stored in the column vector b. 3MATLAB stores matrices in memory such that the order of entries in columns is preserved. Thus A(2) = 2, A(5) = 7, etc. You can drop the geometry and get all the entries of a matrix as a single column vector with the syntax A(:). Sometimes that knowledge or that operation comes in handy. 2 The loop-based way to think about doing this would be to loop over all times, and for each time, loop over all sensor values, multiply each by its weight, and add up those products. Because it would be a pain to type in all the separate statements every time we wanted to do this, let’s write a program. With your favorite text editor, create the following file, and name it sensor.m: M = 5000; N = 500; % Let’s choose M and N reasonably large.A = rand(M,N); b = rand(N,1);C = zeros(5000,1);for i=1:5000for j=1:500C(i) = C(i) + A(i,j) * b(j);endendA few notes about this program: 1. This is a script, a collection of MATLAB commands put into a file with the .m extension; typing the filename is then equivalent to typing in all those commaand separately at the command line. The other type of M-file is a function. Creating those is slightly more complicated; you can read about it by typing help function at the command line. The use of functions is also different from the use of scripts in various subtle ways that I’m not going to get into here, mostly having to do with the scope of variables in memory. 2. In the first line, you’ll see that you can put more than one statement on the same line; semicolons separate statements and suppress echoing, commas separate statements and echo the result. 3. Everything after the % character on a line is treated as a comment and ignored by the interpreter. 4. The rand(p,q) command creates a p × q matrix and fills it with random numbeer between 0 and 1 drawn from a uniform distribution (every value is equally likely). randnhas the same syntax, but draws its random numbers from a Gaussiia distribution with mean 0 and standard deviation 1. Other commands in the same family include zeros and ones, which fill the matrices they create with zeros and ones, respectively. MATLAB’s intuitive that way. Be warned: if you provide only one argument to any of these commands, it won’t give you a vector of that length, as you might expect, but rather a square matrix. 5. for looping is done according to the syntax for = ; ; end. is typically a vector; as the loop executes, takes on each of the values in that vector in turn. So you can have a loop like for i=[6 3 12 -1 7], if you want; if it’s not clear what that does, try it out, having it print the value of i on each iteration. 3 All right, with all that said, we can now run the program. Type sensor at the commaan line. Pay attention to how long it takes to execute. You may have noticed that this situation was constructed such that we can do exacctl the same thing without looping, just by multiplying A and b together (compare ⎡ the definition of matrix multiplication: if C = A � b, then Cij = Aik bkj . Try k C=A*b, and compare its execution time to the looping case. This is why the problem set stresses so much the need to avoid loops; not only does it make your code much more elegant and readable, but MATLAB is full of optimizations to make these kinds of operations on matrices take place much more quickly. 1.3 Plotting and printing One of the great things about MATLAB is the ease and control with which you can use it to make plots. The command help plot will tell you all about it, but I’ll give you a few examples to give you an idea of what you can do: t = 0:0.1:10; % say we take samples every .1 seconds from 0 to 10 secondsy = sin(2*pi*t/5); % make y a sine wave with period 5 secondsplot(y) % the simplest way to plot; notice that the horizontal axis% is not 0 to 10 as we would like, but rather corresponds to% the indices of y, by defaultplot(t,y) % specify both the horizontal and vertical quantities to% plot against one anotherplot(t,y,’go--’) % change the appearance of the plot with an optional% third argument, a string that can specify color% (here, green), point appearance (circles), and% line style (dashed); look in ’help plot’ for more% optionsx = cos(2*PI*x/5);plot(x,y) % parametric plots are equally possibleplot(t,x,’g-.’,t,y,’r*’) % more than one plot can be put on the same% axes, by adding more triplets of arguments% to the plot commandplot(t,x,’k:’) % note that making a new plot by default clears the% old onehold on % after this command, future plots will appear on the same% axes without clearing themplot(t,y,’m-.’)hold off % and now future plots will clear the axes againplot(x,y) semilogx(x+2,y+2) % make the x-axis logarithmic semilogy(x+2,y+2) % or the y-axis loglog(x+2,y+2) % or both axes To print out your masterpiece, MATLAB helpfully provides a print command, which has a longer help file than any other command I can think of. Here are the options you’re most likely to use: 4 print % send current figure to the default printerprint -Pname % send figure to printer called ’name’print -dps file.ps % print to a PostScript file called ’file.ps’print -dpsc file.ps % or color PostScriptprint -depsc2 file.eps % or color encapsulated PostScriptprint -djpeg file.jpeg % or JPEG1.4 How to learn MATLAB It’s actually pretty difficult to go about explaining how to use MATLAB in a methodical way. The best way to learn, and the primary way I learned, is to use its excellent help system, and play around with commands. To get help on the syntax for any command, type help . A lot of commands have intuitive names; if you want to know how to find the mean of a vector, for instance, try help mean, and there it is. If you want a command that performs some function but you can’t think of the name of the command, use lookfor ; for instance, try lookfor variance. Here’s a set of commands that you might find especially useful; I’ll let you read about how to use them. As always, ask any questions you like. while: the other kind of loopif/elseif/else: the other use of conditional statementssize(A) : returns M and N for an M-by-N matrix Arepmat(a,b,c) : tiles matrix A, repeating it b times vertically and ctimes horizontally; try ’repmat(rand(2,2),2,3)’max, min : note that ’m = max(x)’ returns in m the maximum value inthe vector x; ’[m,i]’ = max(x)’ returns the same thing, plus theindex of that value in i!mean, sumfind : _very_ useful; read the help filewho : returns the names of all the variables in memorywhos : like who, but returns more detailed informationclear : clears all variables from memory; ’clear x y z’ clears justvariables x, y, and zwhy : answers any question (try it the next time you get an error!)2 Variance and covariance In class, Sebastian explained that variance is a measure of how much a function or time series fluctuates about its mean. It differs from second moment, ∗x2�, in that the mean is subtracted; thus with x(t) = 10 + sin(t)/10 and y(t) = 2 �sin(t), the second moment of x will be much greater than that of y, but the variance of y will be much greater—which fits with our intuitive notion of what variance should be. Sebastian gave the definition Var[x] �∗x2�−∗x�2 . I find it easier to see how this corresponds to subtracting the mean by starting from a formulation where the mean is 5 subtracted more explicitly: Var(x)= ∗(x −∗x�)2� 2 = ∗(x 2 −2x∗x�+ ∗x�) �� 1 �1 2 1 = xi −2 xi ∗x�+ ∗x�2 NN N ii i � 1 �1 = ∗x 2�−2∗x� xi + ∗x�2 1 NN ii 2 = ∗x 2�−2∗x�+ ∗x�2 2 = ∗x 2�−∗x�. Here we used the definition ∗x� = 1 ⎡ xi, and relied on the fact that ∗x� is a fixed Ni quantity which does not depend on the index i and thus can be taken outside the sum. Covariance is a measure of the extent to which two variables fluctuate together; if both tend to increase at the same time and decrease at the same time, the covariance will be high. You can see that if x and y are each large and positive at the same time as one another, and large and negative at the same time as one another, the term ∗xy� will be large and positive, reflecting the high covariance. However, this term is not enough on its own; for the same reasons as with variance, to avoid getting artificially high or low values for the result, we want to subtract the mean from each variable. Again starting from an expression where we do so explicitly, Cov(x, y) = ∗(x −∗x�)(y −∗y�)� = ∗xy�−∗x∗y��−∗∗x�y�+ ∗∗x�∗y�� 1 � 1 � 1 � = ∗xy�− N xi ∗y�− N ∗x�yi + N ∗x�∗y� i i i 1 � 1 � 1 � = ∗xy�−∗y� N xi −∗x� N yi + ∗x�∗y� N 1 i i i = ∗xy�−∗y�∗x�−∗x�∗y�+ ∗x�∗y� = ∗xy�−∗x�∗y� which is the definition Sebastian gave in class. 3 Linear regression You should be familiar with the idea of solving systems of linear equations. For instaance if you have the equations 3x + 2y =8 and x −4y = −3, and you want to solve for x and y, you can write this in matrix form: ⎤ ⎦⎤⎦⎤⎦ 32 x 8 = 1 −4 y −3 6 � If we define ⎤ ⎦ ⎤ ⎦ ⎤ ⎦ A = 3 1 2 −4 , x = x y , d = 8 −3 , the equation becomes Ax = d, which has solution x = A−1d. In MATLAB, you can find this solution with the command x = inv(A)*d or x = A\d. Geometrically, if you think of x and y as coordinates, we’re finding the location where these two lines meet. Or the same approach can describe a different situation: with the familiar equation for a line y = mx + b, the equations 3=7m + b and 1=3m + b give the matrix equation ⎤ ⎦⎤⎦ ⎤⎦ 71 m 3 = 31 b 1 The interpretation suggested here is that we’re finding the line connecting the two points (7, 3) and (3, 1). If you know much about matrix algebra, you know that things aren’t always this simple. This is a case with a unique solution, one and only one. Sometimes the matrix isn’t invertible, e.g., ⎤ ⎦⎤⎦ ⎤⎦ 12 x 3 = 24 y 6 This is the case where the two lines are identical, or the two points are the same; you don’t have enough information to lock the solution down exactly. In this case the problem is called underconstrained. If your matrix A has fewer rows than columns, it’ll necessarily be underconstrained. If you’ve got two equations and four unknowns, you’re stuck; there’s no unique solution, there’s an infinity of solutions. The opposite situation is when you’re overconstrained. Here you also have no unique solution, but in this case it’s not that you have infinitely many solutions, but rather that you have none; there are too many constraints and they can’t all be satisfied at once. Your lines don’t intersect in a single point, or (and this is clearly the case I’m working up to) your points don’t all lie on a single line. In that case, the best you can do in trying to fit a line to those points is to minimize the total error, for some convenient measure of error. The one that’s always used, and this is the one Sebastian talked about in class, is the total squared error, the sum of the squares of the vertical distances from each point to the line: 1 E =(yi − (mxi + b))2 2 i Minimizing this squared error gives rise to the familiar term ‘least-squares’. For the problem to be overconstrained, A must have more rows than columns.4 Since A isn’t square, you can’t simply take its inverse. But what you can do is multiply both sides of the equation on the left by its transpose: AT Ax = AT d 4More rows than columns doesn’t necessarily mean the problem is overconstrained; it’s a necessary but not sufficient condition. The data set could have a unique solution or be underconstrained. 7 This is called the normal equation. And this is the equation that, for a given A and d, minimizes that squared error. In that sense, it makes the left and right sides of the equation as close as they can get. Now look at this matrix AT A. If A is M × N , then AT A is (N × M )(M × N )= N × N ; so it’s square, which means we can at least try to take an inverse, and in almost all real-world cases an inverse exists. Then the least-squares solution is x =(AT A)−1AT d. −1ATThis (AT A)is called the pseudoinverse. You can find this solution in MATLAB by typing in the expression x = inv(A’*A)*A’*d, by using pinv, or by writing A\d. I think all of these are actually calculated differently, but in real-world cases of real-valued data they should all give the same result. 8