Andrew!
Programming, statistics, music, whatevs.
Thursday, November 6, 2014
This blog is defunct! New location: andrewsturges.com/blog/
Find new content (I hope) at andrewsturges.com. I won't be updating or monitoring this blog anymore.
Friday, March 29, 2013
Exploring the Ruby gem craig
I've decided to take down this post after receiving a cease and desist letter from Perkins Coie LLC, which represents craigslist.org. https://github.com/arsturges/craigslist
Wednesday, March 20, 2013
Exploring the Craigslist Ruby Gem
I've removed this post because I got a cease and desist letter from Perkins Coie LLC regarding a violation of the craigslist terms of use. https://github.com/arsturges/craigslist
Monday, October 15, 2012
Solving an Escape Riddle Using a Monte Carlo Simulation in Python
I recently heard this riddle from a friend at a campfire. I wasn't smart enough to come up with the answer, but when it was explained to me I wanted to know how long it would take the prisoners to be released. In this post I'll describe the riddle and the answer, and show how to use Python to set up an experiment, use a Monte Carlo simulator to run the experiment repeatedly, and briefly discuss some of the results.
But this result depends on the order in which the prison keeper (in our case represented by random.randint) chose the prisoners. Running this code at the command line produces different results each time, usually in the range of 2500--3000 (days). To get a better picture of how long it would take a new group of prisoners to earn their release we'd have to simulate the whole range of possible outputs (or which there are nearly infinitely many) But I don't want to run this manually dozens of times and tabulate the results by hand; I want to make the code run the simulation thousands of times, with random inputs each time, and then aggregate the results. That's the Monte Carlo method.
The last step in the Monte Carlo method is to aggregate and interpret the results. In this case, that means taking the average of all the simulation runs. A simple way to do this with lists is to sum the list :
The Riddle
There are 50 prisoners, and one prison keeper. The prisoners will be placed in separate cells, never able to communicate. Once per day, the prison keeper will randomly select one prisoner and take him into a room with two levers. The levers have two positions each: Up, and down. The prisoner may operate the levers however he wishes, then he must return to his cell. (The levers don't do anything; they just represent a binary state.) The prisoners have one chance to indicate to the prison keeper whether all the prisoners have been in the lever-room. If they so indicate too early, they will never be released. If they're correct (even if they're late), they will be released immediately. The prisoners get one chance to meet prior to imprisonment, to discuss strategy. How can they ensure their release?The Answer
Only one lever is needed. Elect one prisoner to be a 'counter.' When a prisoner enters the lever room and the lever is in the 'down' position and he has never moved the lever, then and only then will he move the lever to the 'up' position. Only the 'counter' prisoner may return the lever to the 'down' position, and each time he does so he increments a counter. When the counter reaches 50, he knows that all prisoners have been into the lever room at least once. He can then inform the prison guard, and all the prisoners will be released.A Solution in Python
I came up with the following model in Python to describe the prison scenario. There are probably better and more optimized way to represent this situation, but here's the first one I came up with.#Iteration 1: The basic scenario import random day_number = 0 prisoner_designee = 1 # chosen by prisoners to be the official counter counter = 0 # as tracked by the prisoner_designee lever_position = 'down' # starting position of the lever prisoners_who_have_operated_lever = [] while counter < 50: prisoner_number = random.randint(1,50) # Select a prisoner if not prisoner_number in prisoners_who_have_operated_lever: if lever_position == 'down': prisoners_who_have_operated_lever.append(prisoner_number) lever_position = 'up' if prisoner_number == prisoner_designee: # If it's the designee... if lever_position == 'up': # ...and the lever is up... lever_position = 'down' # ...put the lever back down... counter += 1 # ...and increment the counter. day_number += 1 print day_numberThe above code, iteration 1, uses
random.randint
to generate a random integer. I used if
statements to describe the logic, and printed out the result. Example:$ python prisoner1.py 2536This means that for this particular run, it would have taken the prisoners (or rather, the prisoners' designated counter prisoner_designee) 2,536 days to count 50 'up' levers, proving that all prisoners had entered the lever room and ensuring their release. That's just under seven years(!).
But this result depends on the order in which the prison keeper (in our case represented by random.randint) chose the prisoners. Running this code at the command line produces different results each time, usually in the range of 2500--3000 (days). To get a better picture of how long it would take a new group of prisoners to earn their release we'd have to simulate the whole range of possible outputs (or which there are nearly infinitely many) But I don't want to run this manually dozens of times and tabulate the results by hand; I want to make the code run the simulation thousands of times, with random inputs each time, and then aggregate the results. That's the Monte Carlo method.
Implementing the Monte Carlo Method
We're going to use a loop to run the experiment, so first, wrap the experiment in a function:#Iteration 2: Wrapping the simulation in a function import random def run_simulation(): day_number = 0 prisoner_designee = 1 # chosen by prisoners to be the official counter counter = 0 # as tracked by the prisoner_designee lever_position = 'down' # starting position of the lever prisoners_who_have_operated_lever = [] while counter < 50: prisoner_number = random.randint(1,50) # Select a prisoner if not prisoner_number in prisoners_who_have_operated_lever: if lever_position == 'down': prisoners_who_have_operated_lever.append(prisoner_number) lever_position = 'up' if prisoner_number == prisoner_designee: # If it's the designee... if lever_position == 'up': # ...and the lever is up... lever_position = 'down' # ...put the lever back down... counter += 1 # ...and increment the counter. day_number += 1 return day_number print run_simulation()Running this code should produce the same type of output (with a different result, of course):
$ python prisoner2.py 2950The next step is to write a loop that will call this function a thousand times:
#Iteration 3: Loop it to run 1000 times import random def run_simulation(): counter = 0 day_number = 0 prisoner_chief = 1 lever_position = 'down' prisoners_who_have_operated_lever = [] while counter < 50: prisoner_number = random.randint(1,50) # Select a prisoner if not prisoner_number in prisoners_who_have_operated_lever: if lever_position == 'down': prisoners_who_have_operated_lever.append(prisoner_number) lever_position = 'up' if prisoner_number == prisoner_chief: if lever_position == 'up': lever_position = 'down' counter += 1 day_number += 1 return day_number simulation_results = [] simulations = 0 while simulations < 1000: simulation_results.append(run_simulation()) simulations += 1 print simulation_resultsIn this iteration we use a
while
loop to run the simulation 1000 times, storing the results in the list called simulation_results. Here's an example of the output:$ python prisoner3.py [2444, 2119, 2818, 2253, 2586, 2543, 2490, 3388, 2034, 2739, 2554, 2585, 2498, 2689, 3180, 2760, 3145, 2698, 2196, 2769, 2400, 2783, 3091, 2258, 2952, 1730, 2974, 2656, 3059, 2000, 3222, 2186, 3114, 2618, .... ](List truncated.)
The last step in the Monte Carlo method is to aggregate and interpret the results. In this case, that means taking the average of all the simulation runs. A simple way to do this with lists is to sum the list :
average = sum(list)/len(list)A more elegant way is to use Numpy's built-in
mean
function, which requires first converting the list
to a Numpy array
:import numpy
array = numpy.array(list)
average = numpy.mean(array)The next iteration of the code will include this small refactorization. Running this code produces the following result:
$ python prisoner3.py 2720.895Now, instead of one experiment, this result represents the average of 1000 experiments.
So what?
So on average, it would take 50 prisoners in this scenario 2721 days to earn their release. Cool. Neat. But this raises so many more questions: Is this number related to the probability of )1/50)*(1/50)? Where is this average value in the range of values yielded from the 1000 results? If we halve the number of prisoners, will the number of days go down by half? What does a histogram of the results look like? The answers, as well as a more in-depth discussion of other statistics, will be in the next post.Wednesday, August 29, 2012
Installing Twitter Bootstrap in Flask 0.9
I have a small web application called 'dispatchtool' written in the Python microframework Flask. I want to use Twitter's Bootstrap CSS framework. In this post I will explain how to do this.
I am running Ubuntu 12.04, using Python 2.7, and Flask 0.9. This tutorial assumes you've got a basic Flask application up and running, and that you can view the development version on your local machine. It also assumes that you're using templates. (This is important!)
To do this, first, create a new directory in your Flask application's root directory called 'static'. My app is called 'dispatchtool', so I create 'static' inside the 'dispatchtool' directory:
Now that we've got our test css file in the right place, we need to tell Flask how to find it. CSS files are always linked to HTML pages in the
Even though that's working, we're not quite done yet. Web application frameworks like Flask have a better way of creating links to static things like CSS files. In Flask, it's a method called
I am running Ubuntu 12.04, using Python 2.7, and Flask 0.9. This tutorial assumes you've got a basic Flask application up and running, and that you can view the development version on your local machine. It also assumes that you're using templates. (This is important!)
Step 1: Get a basic CSS file working.
Before dealing with Bootstrap, we want to ensure sure that the application can correctly find and use a simple test CSS file.To do this, first, create a new directory in your Flask application's root directory called 'static'. My app is called 'dispatchtool', so I create 'static' inside the 'dispatchtool' directory:
$ cd dispatchtool
$ mkdir staticNext, in that new directory, create the test css file. I'll call mine 'test.css'. Add some simple style that will be loud and obvious. It's only purpose will be to ensure that we can get the application to see this file and use it. Here's my 'test.css':
h1 {color: red}
p {color: green}Of course, make sure that your app includes an
<h1>
tag and a <p>
tag somewhere, or else these styles won't be applied, even if the stylesheet is properly set up.Now that we've got our test css file in the right place, we need to tell Flask how to find it. CSS files are always linked to HTML pages in the
<head>
section of each page using the <link … />
tag. In my application, I'll add the following tag:<link type="text/css" rel="stylesheet" href='/static/test.css' />Here's the top of my page, for context:
<!DOCTYPE html>
<html lang="en"> <head> <title>The Demand Response Dispatch Tool</title> <link type="text/css" rel="stylesheet" href='static/test.css' /> </head> <bodyNow open up your page in a browser, and see if the styles were applied. In my case, it suddenly looks like Christmas. If nothing changed, look at the source (ctrl-U) and click on the link you created to the stylesheet. That will show where Flask thinks the CSS file is.
Even though that's working, we're not quite done yet. Web application frameworks like Flask have a better way of creating links to static things like CSS files. In Flask, it's a method called
url_for()
. Methods like this make your application more flexible, allowing you to move around pages without breaking links to your CSS files. So now, replace the value of the 'href' attribute in the <link>
tag with the following line:{{ url_for('static',filename='test.css') }}The double braces tell Flask a) that the contents are to parsed as Python code and b) to render the result directly into the HTML output. Now the whole link tag should look like this:
<link type="text/css" rel="stylesheet" href="{{ url_for('static', filename='test.css') }}" />Now save the file and refresh your browser page. If that still renders your page using your test CSS file, then we're most of the way done. If it doesn't, check the page source, and ensure that the <link> tag is pointing to the right place.
Step 2: Download and "install" Twitter Bootstrap
The next step is simple and short. We need to get the Bootstrap files (download) and put them in the right place (install). Normal people will navigate to the Bootstrap page, click the link, and unzip the files, then put those files (as-is, with no changes to the directory structure) in the 'static' folder that we created earlier. Super elite hackers will do it the following way:$ cd ../static/
$ wget http://twitter.github.com/bootstrap/assets/bootstrap.zip
$ unzip bootstrap.zip
$ rm -r bootstrap.zipWhichever way you do it, you should now have a folder called 'bootstrap' inside your 'static' folder. That's it. You have downloaded and 'installed' Bootstrap.
Step 3: Link to the boostrap.css stylesheet
The main Bootstrap stylesheet is in bootstrap/css/, and it's called 'bootstrap.css'. We want to change our link from 'test.css' to this file, so update your<link>
tag as follows:<link type="text/css" rel="stylesheet" href="{{ url_for('static', filename='bootstrap/css/bootstrap.css') }}" />That's it. If you reload your page, you should see at least some changes reflecting Twitter's style. (When you go to production, you should switch from the general 'bootstrap.css' to the minified 'bootstrap.min.css' to save bandwidth.)
Step 4: Update your HTML to hook into the Bootstrap styles
Now that we've got our Flask app properly linking to the Boostrap css file, all that remains is to use it. That means we have to look at the boostrap.css file, and see what classes and ids we need to assign to our HTML elements. I'll go through two quick examples.Container
The most important is probably to use Bootstrap's "container" class. Simply add the class 'container' the<body>
tag:<body class="container">When you refresh the page, you should see your content contained in a 940 pixel-wide column.
The Hero Unit
The most famous Bootstrap style is probably the "hero unit". To see it, just add the class 'hero-unit' to any block-level element:<div class="hero-unit"><p>Don't abuse this style!</p></div>
Conclusion
We walked through how to set up general CSS file in a Flask app, then how to acquire and link Bootstrap to the same Flask app. The next step is to check out the Bootstrap Documentation, and start choosing the styles you want to use. Thanks for reading!Thursday, July 5, 2012
Basic Unit Testing Tutorial in Python 3.2
I have a Python 3.2 program that's in a working state, but I never wrote any tests for it. Now, I want to go back and provide some code coverage so that as I go back to refactor and add functionality, I'll be more certain I'm not breaking anything. This is a fairly common case--you have existing code, and you want to write some tests for it.
This is now a trivial, but functioning test.
As you can see, when we ran the test suite from inside IDLE, we got the same successful output telling us the test passed (OK!), but it's followed by an ugly error and a traceback. The error is because unittest wants to "exit" but it can't, because the default IDLE behavior is to keep running. There are two simple ways around this error:
In order to call the function in question (
Environment
I'm running Python 3.2 in a Windows 7 environment, using the IDLE Python code editor. I typically run code directly from IDLE by using the F5 key. Because of how IDLE invokes the Python interpreter, some of the code examples below will run differently depending on whether code is invoked from IDLE or from the command line. To enable Python from the Windows command line, it must be added to the PATH. To do this, open a command line window (Start -> Accessories -> Command Line) and typepath <path to your python installation>
. For example, since my python 3.2 installation can be found at C:\Python32, I issue the following path command:>path C:\Python32This adds the Python interpreter to the path so it can be invoked from within your Python project directory. Note that it only changes the path for this one command line session, so if you close and re-open the command line window you'll have to re-issue the path command.
Python's UnitTest Module
The first issue is what testing framework to use. I'm going to use Python's built-in module called UnitTest. It's a one-stop shop for all the testing most pure-Python programs will need. More advanced testing frameworks are available. The two most popular seem to be py.test and nose. I want to explore testing in a more 'raw' state, so I'll save these frameworks until I'm comfortable usingunittest
on its own.Writing Testable Code
It's easiest to write a test for a small, discreet bit of code that accomplishes a well-defined task. Python best practices dictate that code be broken into small functions that perform one task at one level of abstraction. If you follow this advice and write your code using many small functions that each perform one well-defined task, then those functions will be easy to test. If instead you have written large functions that combine many steps, then first consider breaking these functions down into smaller ones.Pick a small function to test
Start small by picking a short, well-understood function that does something simple. I'll use a function that takes one argument, uses anif
statement to choose a "rate", and returns that rate:def choose_renovation_rate(self, years_since_last_renovation): if years_since_last_renovation < 7: rate = 0 elif years_since_last_renovation < 15: rate = 0.01 elif years_since_last_renovation < 25: rate = 0.05 elif years_since_last_renovation < 50: rate = 0.07 else: rate = 0.1 return rate
Create a new file to hold the test code
To get started with the testing itself, create a new file, and for now make sure it's in the same directory as the code you wish to test. I called mineunit_tests.py
. Import the testunit
module, and start a new class that describes what you will be testing. This new class will be a sub-class of unittest.TestCase
. This class will hold several small test cases.import unittest class TestRateFunctions(unittest.TestCase): pass # test cases go here
Write a test function title and description
Next, inside the TestRateFunctions class I just defined, I will define a function and describe in words what it will test. The function names, by convention and to allow the testing machinery to run properly, should start withtest_
. For the moment, I will use a single assertTrue
statement and pass it the condition True
, just to see if everything is working.import unittest class TestRateFunctions(unittest.TestCase): def test_renovation_chooser_should_return_correct_rate(self): # This function should return the correct renovation rate self.assertTrue(True)
This is now a trivial, but functioning test.
Run the test
To actually run the test, we need to add one line to the end of the file:unittest.main()
. With this line in place, we can simply call this test one of two ways.Via the command line
After issuing the PATH commands as discussed under "Environment" above, in the command line terminal navigate to the directory where you created the file and call it:
C:\lighting_floor_space_stock_model>python unit_tests.py . ---------------------------------------------------------------------- Ran 1 test in 0.000s OKThe output consists of
- A period. This represents one test, which in this case was
test_renovation_chooser
. With more tests written, these dots acts as a progress indicator for the test suite. - A long line of hyphens. This is just a visual separator representing the end of the execution of tests.
- A report, which in this case Ran 1 test in 0.000s.
- The phrase OK. We only get this if all tests pass.
Via IDLE
In the IDLE code editor, you can just press F5 to run the code directly. This produces the following output:. ---------------------------------------------------------------------- Ran 1 test in 0.004s OK Traceback (most recent call last): File "C:\lighting_floor_space_stock_model\unit_tests.py", line 8, in <module> unittest.main() File "C:\Python32\lib\unittest\main.py", line 124, in __init__ self.runTests() File "C:\Python32\lib\unittest\main.py", line 272, in runTests sys.exit(not self.result.wasSuccessful()) SystemExit: False
As you can see, when we ran the test suite from inside IDLE, we got the same successful output telling us the test passed (OK!), but it's followed by an ugly error and a traceback. The error is because unittest wants to "exit" but it can't, because the default IDLE behavior is to keep running. There are two simple ways around this error:
- One solution is to just catch the error and pass around it. Simply change the last line to following:
try: unittest.main() except SystemExit: pass
- Alternatively, we could keep the first line as-is (without adding the try/except statements or the second line) and just add an argument to the original line, as follows:
unittest.main(exit = False)
Now Test Your Function
Now that we're sure we've got all the machinery operating correctly (importunittest
, create a subclass, define a test case starting with test_
, make sure it runs properly), we have to write code that actually tests our function.In order to call the function in question (
choose_renovation_rate
) inside the testing function (test_renovation_chooser_should_return_correct_rate
), we need to first make sure its parent class is available. In this case, choose_renovation_rate
is a function inside the class FloorSpace(), so first I need to import the FloorSpare class with from floor_space import *
. Then I can create a new variable called self.rate
that holds the result of the FloorSpace.choose_renovation_rate function. The whole file now appears as follows, using assertTrue
:import unittest from floor_space import FloorSpace class TestRateFunctions(unittest.TestCase): def test_renovation_chooser_should_return_correct_rate(self): # A 2-year old building has a 0% chance of renovation: self.rate = FloorSpace.choose_renovation_rate(self, 2) self.assertTrue(self.rate == 0) unittest.main(exit = False)
A Quick Refactor
Let's make one short change. Since we're testing for equality (using the == operator), we could use the assertEqual function instead of assertTrue. Furthermore, we could a string argument to the end of the function that will get printed if the test fails, allowing us to pass more detailed and useful information to the user. Change the last line of the test function as follows:self.assertEqual(self.rate, 0, "2-year-old building should have a 0% renovation rate.")This is such a simple test method that the message is nearly redunant to just reading the code, but it serves to illustrate the method of adding an information string to the assert method.
Refining the Test
Consider again the function we're testing:def choose_renovation_rate(self, years_since_last_renovation): if years_since_last_renovation < 7: rate = 0 elif years_since_last_renovation < 15: rate = 0.01 elif years_since_last_renovation < 25: rate = 0.05 elif years_since_last_renovation < 50: rate = 0.07 else: rate = 0.1 return rateWe pass it the number of years since it was last renovated (in years), and it returns a renovation rate. That renovation rate is used by other functions to perform various other computations. We could test each of the if conditions, but those numbers might change if we later tune the function to match real-world data. It might be more useful to test, for example, that the function always returns a percentage (i.e. a float between zero and one). So let's test a few of the if conditions just for completeness, and add a new test function to check whether it always returns a percentage between 0 and 1. Note that instead of testing each year we put the assert methods inside a for loop.
import unittest from floor_space import FloorSpace class TestRateFunctions(unittest.TestCase): def test_renovation_chooser_should_return_correct_rate(self): # A 2-year old building has a 0% chance of renovation: self.rate = FloorSpace.choose_renovation_rate(self, 2) self.assertEqual(self.rate, 0, "2-year-old building should have a 0% renovation rate.") # A 13-year old building has a 1% chance of renovation: self.rate = FloorSpace.choose_renovation_rate(self, 13) self.assertEqual(self.rate, 0.01, "13-year-old building should have a 1% renovation rate.") # A 55-year old building has a 10% chance of renovation: self.rate = FloorSpace.choose_renovation_rate(self, 55) self.assertEqual(self.rate, 0.1, "55-year-old building should have a 10% renovation rate.") def test_renovation_chooser_should_return_percentage(self): # The function should always return a rate such that 0 <= rate <= 1 for years_since_renovation in range(100): self.rate = FloorSpace.choose_renovation_rate(self, years_since_renovation) self.assertGreaterEqual(self.rate,0, "renovation rate should be >= 0") self.assertLessEqual(self.rate,1, "renovation rate should be <= 1") unittest.main(exit = False)For a list of assert methods made available to
unittest
, see the documentation. Conclusion
The basic steps to writing unit tests in Python are:- Write testable code by keeping functions short; each function should perform one task
- Create a new file for the test code (e.g.
unit_tests.py
) Import unittest
- Subclass unittest.TestCase (e.g.
class TestRateFunctions(unittest.TestCase)
) - Create new test functions with description names starting with
test_
(e.g.
def test_renovation_chooser_should_return_percentage)
- Make your "work" functions available to the "test" functions by importing where necessary (e.g.
from floor_space import FloorSpace)
- Add string arguments to pass relevant information to the user in the case of a failed test.
- Refine and refactor for readability and usefulness.
Thursday, May 10, 2012
Basic linear algebra with Numpy on Python 3.2
I'm taking Stanford's online course on machine learning through Coursera. The second week has a good overview of linear algebra and matrix operations. The instructor has provided a useful PowerPoint deck in which he explains the basics. I'm going to go through this pdf and implement the linear algebra using NumPy. I have NumPy 1.6 installed with Python 3.2.
Ml SlideMachine Learning: Linear Algebra Reviews Lecture3
First, open an interactive Python shell and load NumPy:
Division is just multiplication by a fraction:
Ml SlideMachine Learning: Linear Algebra Reviews Lecture3
First, open an interactive Python shell and load NumPy:
>>> from numpy import *
Create a Matrix
Let's try just creating the 4x2 matrix he shows in slides 2 and 3.The basic form for creating arrays is to use the array method with parenthesis:a = array()Inside the parenthesis, nest some square brackets, and in those brackets just put comma-separated lists of elements in more square brackets. Each square bracket represents a row. Make sure they all have the same number of columns. So to create the 4x2 matrix from slides 2 and 3 use the following
>>> a = array([[1402,191],[1371,821],[949,1437]])
>>> print(a) [[1402 191] [1371 821] [ 949 1437]]
Matrix Elements
Next he talks about accessing the matrix element-wise, so a11 should access '1402.'>>> a[1,1] 821Instead of returning the first row of the first column, it gave us the second row of the second column. This is because NumPy is using 0-indexing (start at 0) instead of 1-indexing (start at 1). So to get the first row of the first column we index from 0:
>>> a[0,0] 1402
Matrix Addition
Next let's create two 3x2 matrices and add them together. First, instantiate the matrices:>>> b = array([[1,0],[2,5],[3,1]]) >>> c = array([[4,0.5],[2,5],[0,1]]) >>> print(b) [[1 0] [2 5] [3 1]] >>> print(c) [[ 4. 0.5] [ 2. 5. ] [ 0. 1. ]]Add b and c and see what happens:
>>> b+c array([[ 5. , 0.5], [ 4. , 10. ], [ 3. , 2. ]])As you can see, NumPy correctly performed an element-wise addition.
Scalar Multiplication
Next, multiply a scalar by a 3x2 matrix. We'll use matrix 'b' from above:>>> 3 * b array([[ 3, 0], [ 6, 15], [ 9, 3]], dtype=int32)Again, NumPy correctly multiplied each element of the matrix by 3. Note that this can be done in any order (i.e. scalar * matrix = matrix * scalar).
Division is just multiplication by a fraction:
>>> d / 4 array([[ 1. , 0. ], [ 1.5 , 0.75]])
Combination of Operands
Order of operations is important. In this slide the instructor sets up three vectors (3x1) and provides an example in which he multiples by a scalar then adds then subtracts then divides. Let NumPy do it and see what happens:>>> e = array([[1],[4],[2]]) >>> f = array([[0],[0],[5]]) >>> g = array([[3],[0],[2]]) >>> 3 * e + f - g / 3 array([[ 2. ], [ 12. ], [ 10.33333333]])As before, NumPy produces the same answer as the instructor found by doing it by hand.
Matrix-vector multiplication
We can multiply a matrix by a vector as long as the number of columns of the matrix is the same as the number of rows of the vector. In other words, the matrix must be as wide as the vector is long.>>> h = array([[1,3],[4,0],[2,1]]) # 3x2 >>> i = array([[1],[5]]) # 2x1 >>> h * i Traceback (most recent call last): File "<pyshell#54>", line 1, in <module> h * i ValueError: operands could not be broadcast together with shapes (3,2) (2,1)My multiplication operation didn't work, and it helpfully gave me the shapes of the arrays for which multiplication failed. This is because the multiplication operator '*' causes element-wise multiplication. For that to work, the matrices need to be of the same shape (hence the error message shosed us the shapes were different). What we want is the dot product.
The Dot Product
To multiply two matrices using the dot product, use thedot()
method:>>> h = array([[1,3],[4,0],[2,1]]) # 3x2 >>> i = array([[1],[5]]) # 2x1
>>> dot(h,i) array([[16], [ 4], [ 7]])That gives us the correct answer, according to the slides. A longer example:
>>> j = array([[1,2,1,5],[0,3,0,4],[-1,-2,0,0]]) >>> k = array([[1],[3],[2],[1]]) >>> dot(j,k) array([[14], [13], [-7]])This one checks out with the slides as well.
Matrix-Matrix Multiplication
As far as NumPy is concerned, matrix-matrix multiplication is just like matrix-vector multiplication.>>> l = array([[1,3,2],[4,0,1]]) >>> m = array([[1,3],[0,1],[5,2]]) >>> dot(l,m) array([[11, 10], [ 9, 14]])He goes on to give one more example with a pair of 2x2 matrices:
>>> n = array([[1,3],[2,5]]) # 2x2 >>> o = array([[0,1],[3,2]]) # 2x2 >>> dot(n,o) array([[ 9, 7], [15, 12]])Next he provides some context by using the example of home prices. He sets up a 4x2 and 2x3 matrix and multiplies them to quickly come up with price predictions:
>>> p = array([[1,2104],[1,1416],[1,1534],[1,852]]) # 4x2 >>> q = array([[-40,200,-150],[0.25,0.1,0.4]]) # 2x3 >>> dot(p,q) array([[ 486. , 410.4, 691.6], [ 314. , 341.6, 416.4], [ 343.5, 353.4, 463.6], [ 173. , 285.2, 190.8]])
Matrix Multiplication Properties
Show that matrix multiplication is not commutative:>>> A = array([[1,1],[0,0]]) # 2x2 >>> B = array([[0,0],[2,0]]) # 2x2 >>> dot(A,B) array([[2, 0], [0, 0]]) >>> dot(B,A) array([[0, 0], [2, 2]])To test this another way I asserted equality between the two operations and found a neat element-wise comparison:
>>> dot(A,B) == dot(B,A) array([[False, True], [False, False]], dtype=bool)To show the associative property of arrays, create another 2x2 array C and multiply them in different groupings, but in the same order, to show that the result is always the same:
>>> C = array([[1,3],[0,2]]) # 2x2 >>> A = array([[1,1],[0,0]]) # 2x2 >>> B = array([[0,0],[2,0]]) # 2x2 >>> C = array([[1,3],[0,2]]) # 2x2 >>> dot(A,dot(B,C)) array([[2, 6], [0, 0]]) >>> dot(dot(A,B),C) array([[2, 6], [0, 0]])
Identity Matrix
NumPy comes with a built-in function for producing an identity matrix. Just pass it the dimension (numnber of rows or columns) as the argument. Optionally tell it to output elements as integers in order to clean up the output:>>> identity(3) array([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.]]) >>> identity(3, dtype=int) array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])Show that for any matrix A, AI=IA=A. We can use the same A from before:
>>> A = array([[4,2,1],[4,8,3],[1,1,0]]) # 3x3 >>> I = identity(3, dtype=int) >>> dot(A,I) array([[4, 2, 1], [4, 8, 3], [1, 1, 0]]) >>> A array([[4, 2, 1], [4, 8, 3], [1, 1, 0]]) >>> dot(A,I)==dot(I,A) array([[ True, True, True], [ True, True, True], [ True, True, True]], dtype=bool) >>> dot(A,I) == A array([[ True, True, True], [ True, True, True], [ True, True, True]], dtype=bool)
Inverse and Transpose
Show that if A is an mxm matrix, and if it has an inverse, then A(A-1) = (A-1)A = I. To do this, we can use the same 3x3 matrix A from above:>>> A = array([[4,2,1],[4,8,3],[1,1,0]]) # 3x3
>>> inv(A) Traceback (most recent call last): File "<pyshell#112>", line 1, in <module> inv(A) NameError: name 'inv' is not definedWe got this error because though we loaded NumPy, we need to also load the special linear algebra library:
>>> from numpy.linalg import *Then we can try the inversion again:
>>> inv(A) array([[ 0.3, -0.1, 0.2], [-0.3, 0.1, 0.8], [ 0.4, 0.2, -2.4]])Now show that that A(A-1) = (A-1)A = I:
>>> dot(A, inv(A)) array([[ 1.00000000e+00, -2.77555756e-17, 0.00000000e+00], [ -2.22044605e-16, 1.00000000e+00, 0.00000000e+00], [ -5.55111512e-17, -1.38777878e-17, 1.00000000e+00]]) >>> dot(inv(A), A) array([[ 1.00000000e+00, 0.00000000e+00, -5.55111512e-17], [ -2.22044605e-16, 1.00000000e+00, -5.55111512e-17], [ 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])Note that this was meant to return the identity matrix, but that the float operations returned not zeros but numbers very close to zero. Note also that a matrix of all zeros has no inverse:
>>> C = array([[0,0],[0,0]]) >>> C array([[0, 0], [0, 0]]) >>> inv(C) Traceback (most recent call last): File "<pyshell#121>", line 1, in <module> inv(C) File "C:\Python32\lib\site-packages\numpy\linalg\linalg.py", line 445, in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype))) File "C:\Python32\lib\site-packages\numpy\linalg\linalg.py", line 328, in solve raise LinAlgError('Singular matrix') numpy.linalg.linalg.LinAlgError: Singular matrixThe error message agrees with the machine learning instructor, who calls these special matrices "singular" or "degenerate."
Matrix Transpose
Lastly, show some matrix transposition, whereby the rows are flipped element-wise:>>> A = array([[1,2,0],[3,5,9]]) >>> A.transpose() array([[1, 3], [2, 5], [0, 9]])If we name the transposed array 'B', then we can show that the ijth element of A is the jith element of B:
>>> B = A.transpose() >>> A[0,2] == B[2,0] True
>>> A[1,2] == B[2,1] True
Conclusion
We used NumPy and NumPy's librarylinalg
to go through the linear algebra review slides from Coursera's Machine Learning course.
Subscribe to:
Posts (Atom)