Computation Frameworks for the Social Sciences (aka I’m teaching a class)

Posted on by
Reply

This spring I am teaching my first course!  It is a pretty small seminar, 8-10 graduate students, but they all seem excited about the material.  It is a “Programming for Political Scientists” course that will use Python to both teach people how to write good software as well as show them how people are using software in the discipline currently.  I hope to spend the first half of the course covering basic software engineering and computer science concepts before moving on to some specific applications.  Hopefully, by the end of the course all of the students will have built something that they can use to further their research agendas (e.g. a web scraper to supplement a data set).

The class website is up at http://joshcutler.github.com/PS398/ where you can find my schedule outline as well as the homeworks.  One thing that has worked well so far was requiring that everyone work through Zed Shaw‘s Learn Python the Hard Way before class.  This got everyone (note that most had not done any serious programming before) on the same page and ready to start tackling some more advanced ideas.

For those of you that have experience in the space, I would welcome any feedback on the syllabus (or anything else).  As the course moves along I will try to post about any changes that I make, findings, etc. for those interested in teaching a similar course elsewhere.

Line numbers on embedded Gists

Posted on by
1

For all of you bloggers out there that like embedding gists but are frustrated by the lack of line numbers, I found a nice CSS solution.  After a little googling, I found this solution by potch.  It looks as though the structure of an embedded gist has changed slightly so I modified the css to look like the following:

.gist .highlight {
    border-left: 3ex solid #eee;
    position: relative;
}
 
.gist .highlight pre {
    counter-reset: linenumbers;
}
 
.gist .highlight pre div:before {
    color: #aaa;
    content: counter(linenumbers);
    counter-increment: linenumbers;
    left: -3ex;
    position: absolute;
    text-align: right;
    width: 2.5ex;
}

Include this somewhere in your stylesheets and viola, you get nice line numbers for all of your embedded gists (in modern browsers at least).

Rook

Posted on by
Reply

As I said in my last post, I am working on setting up a statistical web service using R.  I decided to use the Rook library to do so and wanted to give a brief overview of how Rook works for others who might be interested.

The best way to learn is often to just dig right in and go line by line, so here is a simple rook application:

require('Rook')

library(Rook)
library(rjson)

rook = Rhttpd$new()
rook$add(
  name ="summarize",
  app  = function(env) {
    req = Rook::Request$new(env)
    
    numbers = as.numeric(unlist(strsplit(req$params()$numbers, ",")))
    
    results = list()
    results$mean = mean(numbers)
    results$sd = sd(numbers)

    res = Rook::Response$new()
    res$write(toJSON(results))
    res$finish()
  }
)

rook$browse("summarize")

view raw rookexample.R This Gist brought to you by GitHub.

Line 6 creates a new rook webserver object.  This is what will respond to HTTP requests from the browser.

Line 7 adds an “application” to the webserver.  When you add an application to a Rook server you need to name a route and specify what should happen when a user requests that action.  The route is specified by line 8 and tells the Rook server that when a user requests “summarize” it should execute the code specified in lines 9-21.  For some reason Rook prefaces all routes with the word “custom”, so the url to access the route we specified would be http://server:port/custom/summarize

The interesting part of the application occurs in the function that is assigned to app on line 9.  Rook wraps the parameters of the HTTP request in some nice accessors which can be used as seen in line 10.  While the docs specify all of the information available, the important part of the HTTP request for our app is the params() method.  This returns the union of any variables passed via query string and POST.  For the rails developers it is exactly the Rails params hash.

Line 12 parses the input.  This application is expecting an array of numbers separated by commas assigned to the numbers parameter.  Note that if this isn’t specified the app will break in its current state.

After parsing the input parameters we compute some simple summary statistics on this array of numbers and store them in the list results (14-16).

Rook uses the Rook::Response object to fashion a proper HTTP response.  After instantiating the response on line 18, we call the write() method on the response object.  This sends whatever string it is passed to the output stream which is returned to the requestor (in our case the web browser).  Note that I am returning the results in JSON format and should probably set HTTP headers as well if I was going to deploy this to production.  Line 20 flushes the output stream and returns it to the requestor, call this when you are done constructing your response.

Finally, to start our server and see the thing in action we can use the browse() method (line 24).  We are specifying the action that we want to browse to as the parameter.  This should pop open a web browser pointing to your Rook application–and you should get an error:

Error in strsplit(req$params()$numbers, ",") : non-character argument

This is because we didn’t validate the input and our application is expecting an array of numbers.  So, lets pass them through the query string.  Append the following to the URL in your browser: “?numbers=1,2,3,4,5″ and refresh.  No you should get the following:

{"mean":3,"sd":1.58113883008419}

And you have a functioning Rook server!

If you want to add other functionality you can just make addition add() calls on the rook object and the new applications will be added.  If you want to change the functionality of an existing application make sure that you remove (rook$remove(‘summarize’)) it and then re-add it (or just restart everything).

Rook and R Webservices

Posted on by
1

Recently I have been working on setting up a webservice that does some non trivial statistical work.  Normally, my go to when building web services is Ruby/Rails due to ease of use, then I offload anything computationally intensive to something more optimized (e.g. a C or Java application on the same box).  In this case however, partly because of my co-author’s skill, partly discipline norms, and a whole lot of R being awesome for this sort of thing, the statistical work is going to be done in R.

While it would have been possible to still build the webservice wrapper in Ruby and then either use one of the existing Ruby wrappers for R (or even spinning up an R process on its own), I wanted to see if I could build the whole thing in R.  As is almost always the case, I was not the first person to think of this, and most of the hard work has already been done.

Rook is an R package on CRAN written by Jeffrey Horner.  For those of you familiar with Ruby and Rack, Rook is very similar.  It provides an interface, using R’s built in rApache webserver, to handle http requests, handle routing, etc.  The more I read about it, the more I was convinced that it was a great solution to the problem that I was working on.  Now if only I get it to run on Heroku….

Well, running Rook on Heroku was surprisingly simple thanks to Noah Lorang’s example which you can find here.

So, how do you get started?  Almost everything that you’ll need is in the README in Noah’s repository, but there are a couple of tricky things to note.

First make sure that you have a Heroku account.  It is free to sign up and you get one free full-time process per project (one single web server in our case).  There are numerous resources (including their excellent help files) to get you through this part.

Next you can either walk through the instructions in Noah’s example (which I ended up doing), or you can do the much easier thing by cloning his repository.  If you do this, then you should be able to deploy it directly.

After you get a running instance of Rook, you will want to write some R code.  To run your own custom R script, replace the “/app/demo.R” in rackup file (config.ru) with the path of your script.  Otherwise, you can just put your code in demo.R.

Because Heroku is a read-only file system, you will need to include any R packages in your source tree so that they are “installed” when you deploy.  Initially you will just have R and Rook installed (if you cloned the existing project).  Because some packages require native compilation, you really should do that compilation on one of heroku’s servers.  In order to do this, you need to ssh into your app server:

heroku run bash

Once you are in the bash shell on your app server you can load R.  From there install any packages that will be dependencies for you project.  When you exit the R shell, do not log out of the ssh session.  If you do the app will reset and you will need to start over (remember it is read only so the file system changes persist only as long as your session).  You then need to figure out how you want to get these changes off of your heroku instance.  First zip them up (I zipped up the whole bin directoy):

tar -cvzf mybin.tar.gz ./bin

Then you can either scp it off of the machine (as Noah suggests):

scp mybin.tar.gz me@myserver.com:~/myproject/bin.tar.gz

Or, if you do not have access to a destination that you can scp (heroku does not have ftp installed), you can do the roundabout method of setting up github as a remote, checking in the tarball, pushing it to github and then cloning that repository locally.  Once you have the tarball on your machine, just untar it into your repositories bin directory, checkin, and deploy.  You will now have access to those packages in R.

I’ll writeup how to actually use Rook in my next post.

 

Naive Bayes with Laplacean Smoothing

Posted on by
Reply

In aiclass.com, we just covered Naive Bayesian Classifiers, and it couldn’t have been more perfectly timed.  Prior to that lecture series, one of the projects that I am working on required that I build a classifier for a large body of data that was getting funneled into the system.  I spent quite a bit of time searching for the best way to do this, hoping that there would be a rubygem that could save me some effort, but much to my chagrin, nothing quite fit the bill–so I started in on building my own.

The basic idea behind a Naive Bayes classifier is that we have some set of documents that have been categorized (into n categories) and want to use this information about our existing labeled documents to predict the category of new, not yet labeled, documents.  It is a pretty direct use of Bayes rule and is probably best understood through an example.

Say you have 5 documents:

  • {subject: ‘Must read!’, text: ‘Get Viagra cheap!’, label: ‘spam’}
  • {subject: ‘Gotta see this’, text: ‘Viagra.  You can get it at cut rates’, label: ‘spam’}
  • {subject: ‘Call me tomorrow’, text: ‘We need to talk about scheduling.  Call me.’, label: ‘not spam’}
  • {subject: ‘That was hilarious’, text: ‘Just saw that link you sent me’, label: ‘not spam’}
  • {subject: ‘dinner at 7′, text: ‘I got us a reservation tomorrow at 7′, label: ‘not spam’}
We have 2 spam message, and 3 real messages.  Each of these messages has a subject and some text that we can use to train our classifier.
Given a new message:
  • {subject: ‘See it to believe it’, text: ‘Best rates you’ll see’, label: ?}
What is the probability that it is a spam message?  Using Bayes rule we can compute it in the following way:
P(spam|subject,text)=\frac{P(subject,text|spam)*P(spam)}{P(subject,text|spam)*P(spam)+P(subject,text|notspam)*P(notspam)}
All of these values can be computed by inspection of the previous documents:
P(spam)=\frac{2}{5}
P(notspam)=\frac{3}{5}
Note that in the case of Naive Bayes you assume independence of your variables (which is probably not true given that the english language is structured).
P(subject,text|spam)=P(subject|spam)*P(subject|spam)
P(subject|spam)=\prod_{word \in subject}P(word|spam)
So for example, in the document we want to classify:
P('see' \in subject|spam) = \frac{1}{5}
You will note that the document to classify has some words that are not in any of the existing classified documents (e.g. ‘believe’).  This will give those conditional probabilities a value of 0, thus making the numerator 0 even though there is definitely a greater than 0 chance of this item being classified as spam.
The solution to this problem is known as Laplacean Smoothing.  In order to perform smoothing you pick some parameter k.  In our case we can set k=1.  This smoothing parameter is added to all probabilities as they are calculated and a normalizing constant is added to the denominator to make it a valid probability.
Thus, with a smoother of size 1:
P('believe' \in subject|spam)=\frac{0 + 1}{5+5*1}
Where does the 5*1 come from in the denominator?  Well we have a smoothing factor of 1 and we have 5 different known values for words in the subject, thus in order to make the known values a true probability distribution we need to add that to the denominator (so it sums to 1).
Like I said, that math here is pretty straightforward if you can buy the assumptions.  And, even if you can’t it seems to work pretty well.
So, how did I end up using it in my app?  I built a pretty simple gem to do classification called Classyfier.  It was based loosely on Bayes Motel but I cleaned up and reorganized some things (as well as added smoothing). I anticipate that I will be adding more features to this package as my need for more sophisticated classification grows.  For more info on how to use the gem see the example below or just checkout the test file.
require 'classyfier'

@classyfier = Classyfier::NaiveBayes::NaiveBayesClassifier.new
@classyfier.train({:subject => 'Must read!', :text => 'Get Viagra cheap!'}, :spam)
@classyfier.train({:subject => 'Gotta see this', :text => 'Viagra.  You can get it at cut rates'}, :spam)
@classyfier.train({:subject => 'Call me tomorrow', :text => 'We need to talk about scheduling.  Call me.'}, :not_spam)
@classyfier.train({:subject => 'That was hilarious', :text => 'Just saw that link you sent me'}, :not_spam)
@classyfier.train({:subject => 'dinner at 7', :text => 'I got us a reservation tomorrow at 7'}, :not_spam)
        
@scores = @classyfier.classify({:subject => 'See it to believe it', :text => 'Best rates you\'ll see'})

view raw gistfile1.rb This Gist brought to you by GitHub.

 

Exponential growth and resource consumption

Posted on by
Reply

In August, while one the road moving to Duke, I heard a great interview on NPR with David Suzuki.  During this interview he talked about the relationship between anything that grows exponentially and the resources that it consumes.  In particular he focused on how we are likely to misperceive the remaining level of resources because humans are bad at understanding exponential growth (see anything by Kurzweil).

I wanted to retell his analogy because it has really stuck with me and then provide a visualization that I found helpful.

Imagine 2 bacteria in a petri dish.  The petri dish is many, many times bigger than they are.  To keep the math simple lets say that it is 1 billion times larger.  The bacteria in the dish look around and realize that they have a bounty of agar on which to live and thus anticipate much prosperity.  So the bacteria, being bacteria, replicate.  After 1 unit of time (lets say they replicate every minute), there are now 4 bacteria.  These 4 bacteria look around and see that they have this vast petri dish at their disposal and thus keep replicating.  Since replicating is fun, this goes on for a while.  Eventually after 20 turns of doubling in size there are 1,048,576 (2^{20} ) bacteria.  They realize that they have used up just more than 0.1% of their resources and there is much rejoicing and still more replicating.

By the time that 29 minutes have passed there are now 536,870,912 bacteria.  They have used over half of their resources, but they look around and see that they still have half of them left and it took them 29 turns to get here, so they have to have a little time to figure out how to get out of this mess.  On the 30th turn, the population doubles again and now there are 1,073,741,824 and they are out of room.  When a population grows by doubling in size, it by definition will go from 50% resource utilization to 100% in one time period.  That was both obvious (mathematically) and shocking to me.  How fast is the human population growing?  How close are we to running out of resources on earth?

Well, if there was any doubt about whether human population growth was exponential here is a graph of estimated growth over the past 12,000 years (courtesy of wikipedia):

Human Growth Chart

And the CIA factbook estimates the current population growth to be about %1.092 year over year.

Here is a graph that just visualizes the thought experiment proposed by Suzuki:

The point of this graph is just to show how quickly the resource level can plunge to zero and how as an inductive species it is easy to fall into the trap of not understanding the rate at which exponential growth can take off by using historical data to improperly predict the future.

Note: For the graph I just picked a time period of 1000 and then decided to drive it to 0 at that point.  Given that there is a finite amount of resources you will get the same graph form, one way or the other after fixing the resource size.  Below is a gist in R to play around with and see how the curve forms change under different growth rates and how the warning time (as a percentage of the species’ history) shrinks as the size of the resource goes up.

rm(list = ls()) #Clear the workspace
library(ggplot2)

time = seq(1:1000)
starting_population = 2
reproduction_rate = 1.0192

population = cbind(time, starting_population*reproduction_rate^time)

starting_resources = population[time==length(time)-1][2]
resources = cbind(time, starting_resources - population[,2])

data = as.data.frame(rbind(population, resources))
# data set indicator 
data$ds <- factor(rep(1:2, c(length(time), length(time))))
data$Legend <- factor(rep(c("Population", "Resources"), c(length(time), length(time))))

plot1 = ggplot(data, aes(x=time, y=V2, group=ds, color=Legend)) + 
  geom_line(factor=data$ds) +
  scale_y_continuous('Count', limits=c(0, starting_resources)) + 
  scale_x_continuous("Time")
print(plot1)
view raw gistfile1.r This Gist brought to you by GitHub.

AI-Class.com

Posted on by
Reply

When I found out about ai-class.com, this fall I was really excited. It is an online Intro to AI course taught by Sebastian Thrun and Peter Norvig out of Stanford University.  It is free, they are smart and I hadn’t thought about AI problems since undergrad.  Well, it has finally started–and so far so good.

There were a few technical glitches as things got rolling, but that was mostly due to the insanely high demand for the course.  At one point, according to their twitter account, they were getting over 7000 web requests per second.  They now claim to have over 160,000 students registered for the course.  As of today everything seems to be running very smoothly.

When I first heard the numbers regarding how big the course was I wasn’t sure how they were going to administer homework at that scale.  Well, today I turned in my first homework and it is rather cleverly done.  All of the lectures are short video clips (1-6 minutes so far), and at end of each of them they pose some sort of quiz question.  They way that they are presented is pretty clever however, in that they draw out the question in the video via pen and paper and then superimpose an html form on top of it so that you can submit your responses.  For example, here is a multiple choice quiz question:

This ends up working surprisingly well.  So, for the homework assignments they do the same thing.  One of the presenters draws out a question (a maze, or graph, etc) and then asks some questions about it, the form gets superimposed and you submit your answers.  Thus far it has been really great and I highly recommend that anyone who is interested try it out next time the courses are offered.  They also have a Machine Learning course being taught by Andrew Ng as well as a Database course by Jennifer Widom all being done in the same format.

Keeping out the Riff Raff

Posted on by
Reply

I was just wrapping up my MPSA application for a paper that I am working on with Jacob Montgomery and found out that first year graduate students can’t apply to present papers!  It ended up being OK because Jacob was able to submit it for us, but I just wanted to give anyone else out there a heads up.  The deadline is tomorrow and you will need someone more senior to do it if you are planning on submitting a paper.

Getting into Git

Posted on by
Reply

This monday, I was asked by Mike Ward to give a talk to his ICEWS group at Duke about the wonders of source control and how to get started using git.  I was excited to do it because I am a big fan of git and use it on almost everything that I create myself (yes I use it on my CV…).  While the opportunity to introduce it to others was great, I went into it a little bit of trepidation because in the past when I have talked about git or other DVC systems, my audience is at least sold on the idea of version control… because they are software developers.  That was not the case for this audience, but I think I may have converted them nonetheless.

Given the audience, I chose to skip over most of the technical reasons why git is nice.  If you don’t have experience with other forms of source control then it is tough to get sold on the particulars of git (distributed, easy branching, etc.).  This group did however have a code base that multiple people were working on and currently sharing through dropbox. So, I wanted to at least sell them on how Version Control and a distributed host (like github) would make their lives better.

We started by walking through the steps of:

  • Initiating a repository
  • Committing a change
  • Pushing it to a remote host
  • Cloning that repository
  • Pulling changes

At this point one of the audience members asked the obvious question about “what happens when people commit conflicting things?”  This had been a major problem when there was only one instance of the codebase and Dropbox was syncing it.  This was a great question because it segued perfectly into resolving merge conflicts.  We did one of these by hand to show how git handles that and then spent some time looking at github and its nice presentation of commits, branches, etc.  This really seemed to resonate with the audience.

We did not get to branching during the talk, so I wanted to follow up here for those that were interested.  Branching is a way of making a parallel version of the codebase.  If for example, you wish to make some major changes to a library but need to keep the main branch of the build working because others will be editing it as well, you create a branch.  This parallel version of the codebase can then be edited by you.  You can commit changes and even push them to remotes because this is all happening in parallel to the master branch.  When you are done with your work you can merge your changes back into your master branch.  This takes the commits that you made on your branch and applies them to master (very much the same process as a git pull command).  This probably best understood visually so checkout the great explanation on learn.github.com for more.

Any feedback on how to make this presentation more relevant to Political Science would be appreciated.  My goal was to really focus on how it could help in the immediate term and not get bogged down in the technical reasons for why git is good (of which there ware many).  I was happy to see that by the end of the talk a few of the audience members had already setup repositories for both their R codebases as well as their papers that they were doing in LaTex.  Here is my slidedeck, nothing groundbreaking here, but a reference nonetheless:

Intro to git

View more presentations from josh_cutler.