Big data’s all the rage, but sometimes a couple thousand human-generated labels can be pretty effective as well. And since I’ve been using Amazon’s Mechanical Turk system a lot recently, I figured I’d share some of the things I’ve learned.

What is MTurk?

Mechanical Turk is a crowdsourcing system developed by Amazon that connects you to a relatively cheap source of human labor on the fly.

For example, suppose you have 10,000 websites that you want to classify as spam or not. To get these classifications, you (the Requester):

1. Create a CSV file containing the links and any other information.
2. Log onto MTurk and create a HIT (Human Intelligence Task) describing the job (possibly by using Amazon’s WYSIWYG editor or writing your own HTML, which can refer to columns in your CSV). [There’s also an MTurk API, if you don’t want to use the terrible UI.]
3. Within hours of starting the task, your judgments will be completed by Turkers around the world for pennies each.

More Example Tasks

So what can you use MTurk for? Here are three of my favorite uses:

• A Sequence of Lines Consecutively Traced by Five Hundred Individuals

• The Sheep Market: asking Turkers to draw sheep

• Blurry Text Transcription (Seriously! How is this possible?!)

And here are some more practical tasks, from HITs running right now:

• Categorize the sentiment of a tweet towards Panera Bread

• Copy text from a business card

• Judge entity relatedness

Increasing the quality of your judgments

So what will the quality of your judgments look like?

If you don’t do anything special, then your output will contain a lot of garbage. I’ve thrown out entire tasks because of scammers who spend less than 5 seconds on each judgment (Amazon records the time each worker spends) and submit random clicks as output (e.g., labeling Nike as a food category).

Luckily, Amazon provides a few worker filters:

• You can require that only Turkers who have received at least (say) 99% approval rate on at least 10,000 judgments in the past are allowed to work on your judgment. (If you see bad judgments from a worker, you can reject them and get your money back.)
• About a year ago, Amazon launched a “categorization masters” and “photo masters” program, which allows only masters to work on your HITs. According to a chat with a member of the MTurk team, Amazon assigns these master badges by creating special tasks (anonymously, and for which Amazon already knows the answer) and measuring the quality of each worker’s response to these tasks.
• You can also create a custom filter and handpick who gets allowed to work for you, or set up a qualification test that workers are required to take before working on your tasks.

I’ve used different combinations of the first two filters, and gotten excellent results – compared to in-house judges I’ve worked with in person and paid \$20-30 an hour, the judgments on Mechanical Turk have been just as good and sometimes even better. (I often ask my judges to explain their judgments, which makes it easy to detect high quality workers.) For example, here are some typical response I’ve received when asking judges to determine which of two products a given Twitter user might be more interested in:

The user is a female obsessed with Twilight Movies and Rob Pattinson. She tweets and follows both subjects. Movie tickets would be interesting to her.

He doesn’t seem to play video games, and he doesn’t seem technical enough to care about running Windows on a Mac. Neither of these products are a good fit for him.

In fact, I’ll frequently also get emails from Turkers giving me suggestions on how to improve my tasks or asking how they can do them better. (Amazon allows workers to email you. The only way for the requester to initiate a conversation, though, is by paying the worker a small bonus for excellent work, and including a message with the bonus.) Here are excerpts from some emails I’ve received:

I just wanted to check in to be sure that once I figured things out that I was doing your hits the way you intended them to be done. I want to be sure that you are getting the data that you need from the work. Please do not hesitate to let me know if there is anything that I can do to improve the way I am working your HITs. This is my full time job while I stay at home with my kids, so I like to check with the requesters to be sure that I am putting out the work that they are looking for. Any suggestion is welcome.

Frankly, lingerie, makeup, and feminine hygiene are the only male-exclusionary topics I can think of, and it feels knee-jerk sexist to mark any sports-related site for men. That said, should I hew more closely to gender stereotypes or be politically correct? (from a HIT where I was gathering gender classification data)

I do think a few more categories are needed but keeping the number down overall is good - 50 or 60 to choose from can be overwhelming and not worth the time. I may have mentioned I never used the Photography one (and I did a lot of those) so that is a good candidate for elimination.

That said, despite the approval rate filters and masters badges, I do occasionally get a couple scammers in the mix (or even just judges who don’t produce as excellent work). So one suggestion is to run an initial task with these filters applied, find the workers with the best quality, and from then on use a custom pool containing these Turkers alone.

How much to pay

So how much should you pay your workers?

New Turkers and Turkers who don’t meet the strict filters can be paid less, but most of my high-quality workers expect to make about \$8-14 an hour. (You can only specify how much you pay per judgment, but Amazon will tell you how long each item ends up taking on average.) For example, here’s what several Turkers said what I asked them directly how much they make:

Most of the work I do is either writing or editing. When editing work is available, I make \$15-20 per hour. I’m a slower writer than an editor, so I average \$10-12 per hour with writing. I also judge sentiments of messages and average about \$8 per hour with that type of work. I would like to average a minimum of no less than \$8 per hour.

A big factor in deciding to do a task or not comes from the time investment involved. The two big time sinks are either googling/searching/having to go to another site, or having to write something as part of your reply. If I remember correctly, a) your tasks did require looking at another page but either the link was right there OR, better yet, you had that page embedded in the HIT itself so clicking out of the window wasn’t necessary (turkers get very excited about this), and b) the quality of the pay rate was such that it easily outweighed the time it took to leave an explanatory comment.

For me at least, those things can’t be underestimated. Sure, your tasks may be a little time-consuming, but I figure a good task is one I can make 10 to 12 cents a minute on. Your task might take longer but I’m definitely still coming out ahead.

From my own experience, I work hardest and best for a requester that pays well and doesn’t reject (or at least seems to have a reason for a rejection when it happens). If a requester is going to accept the majority of my work, I as a worker feel that obligates me to provide them with the best quality possible. Similarly, although I’m conscientious with all tasks, I’m especially so with a high-paying one: it would be easy to take advantage of a high-pay, low-reject requester - which would ultimately lead them to either lower the pay or change the acceptance criteria. I don’t want that!! That’s the kind of requester I want around. I’m grateful for high pay and fair policies and that kind of requester gets an above-and-beyond effort from me.

For the pay, I have worked on master’s hits that have ranged from \$6-\$16 per hour. Averaging them out works out to around \$9, which isn’t a bad wage. I have two requesters that I work for that don’t use the master qualification but instead have closed qualifications that they’ve assigned to their best workers. Those tasks pay between \$12 and \$15 per hour, so no matter what I’m working on I will stop what I’m doing to work on them. The best paying hits are always done very quickly, so most of the time if you check out mturk and look at the tasks available you won’t get a very good idea of average pay because the terrible paying hits will sit on the board until they expire.

Obviously, this is self-reported, so there’s a strong possibility that the Turkers are artificially inflating their numbers. But this does match what I’ve been told by a manager on the MTurk team, as well as what Turkers self-report on TurkerNation.

A good suggestion regarding pay is to start at the lower end of the scale, around \$6-8 per hour, and increase that until you get both the quality and speed you want.

Other design tips

Interestingly, according to what Turkers (see the excerpts above) and my Amazon contact say, as well as other research I’ve seen (e.g., this paper), pay is not at the absolute forefront of Turkers’ minds when they decide what to work on. Instead, they focus more on requesters they’ve already established a good relationship with, HITs with many items (so they can quickly settle into a rhythm), HITs they know they’ll be paid for (so they’re not worried about rejections), and HITs that they generally enjoy doing more.

So here are a couple suggestions:

• If your task is hard and there’s no clearly correct answer, even good Turkers might be worried that you’ll reject their judgments (and so they might skip over your HIT). So make it clear in your instructions that you won’t reject any judgments, or that you won’t reject any judgments with an honest effort.
• Make your instructions collapsible, or link to them in a separate site. Scrolling is kind of annoying on Mechanical Turk (I know – I’ve tried working on HITs myself), so you should minimize the amount workers have to scroll. Ideally, everything fits on a single screen. Plus, the less workers have to scroll, the faster your HITs will get done. For example, here are excerpts from emails I received from two different Turkers when I first started out:

I have a suggestion that would really make things go a little quicker. Is there anyway you could script the twitter link to automatically open in a new tab? It amazes me how much it can slow you down to have to right click and open it manually in another tab, and when you forget, you have to take a few more steps to get back to where you were.

It would be amazing if the Twitter account could be on the same page instead of having to click to get to another screen - the work would go *exponentially* faster! Overall, I’m enjoying them - and I’m not the only one. Despite your stringent requirements these are disappearing pretty quickly.

• Introduce yourself on TurkerNation, a forum where Turkers and Requesters go to talk about Mechanical Turk. This helps establish your reputation as a good requester who listens to feedback, which will make good Turkers want to work for you. (More on this below.)
• Approve judgments quickly: Turkers want money now instead of money later. For example, one worker told me:

Quick approval is important, too. Watching that money pile up is a serious motivator; I’ll sometimes choose a lower-paying task that approves in close to real time over a higher-paying one that won’t pay out for several days.

When using my trusted set of workers, I let Amazon auto-approve all judgments within a couple hours.

Reputation

Reputation is pretty important. Turkers love requesters who take the time to respond to emails and incorporate suggestions. Excerpts from emails I’ve received:

I LOVE it when requesters care enough to ask the opinion of us lowly turkers and am more than willing to take a few minutes to help them with anything. I look forward to seeing what you cook up!

Thanks for taking the time to try to make your hits better in both pay and design. It’s great to see a requester that actually cares, when most don’t. If you have any other questions for me, feel free to ask. I hope to work for you again soon.

From my own experience, I work hardest and best for a requester that pays well and doesn’t reject (or at least seems to have a reason for a rejection when it happens). If a requester is going to accept the majority of my work, I as a worker feel that obligates me to provide them with the best quality possible. Similarly, although I’m conscientious with all tasks, I’m especially so with a high-paying one: it would be easy to take advantage of a high-pay, low-reject requester - which would ultimately lead them to either lower the pay or change the acceptance criteria. I don’t want that!! That’s the kind of requester I want around. I’m grateful for high pay and fair policies and that kind of requester gets an above-and-beyond effort from me.

I’ve gotten great suggestions from a lot of Turkers (sometimes, when launching a new type of experiment, I’ll do a quick trial run in order to get some fast feedback before spending more time on the HIT design), and I suspect it’s partly because I’ve taken the time to connect with my workers.

So, as suggested above, one way of quickly garnering some goodwill when you’re first getting started is to make a post introducing yourself on TurkerNation. (There’s a sub-forum devoted to this exact purpose, in fact.)

This is useful because workers will often start new threads recommending particular requesters and encouraging other Turkers to work for them. In the amusing thread praising me, for example, one worker mentioned that she’d been hesitant to work on my HITs until she saw the post confirming I was a good requester.

Also, many Turkers mention that they always refer to Turkopticon, a Firefox extension that displays ratings of requesters by other Turkers, before accepting work from a requester they haven’t worked for before.

This is what TurkOpticon looks like:

Here are some comments about TurkOpticon on TurkerNation:

I think that it is well worth taking the time to check reputation of requesters via TurkOpticon and/or in this forum. Checking first substantially minimizes your risk of rejection, of being blocked, and of being paid sub-human wages.

Blindly doing hits for requesters that were never heard of before got me with a pretty bad approval rate when I first started turking. After that, I rigorously inspect every requester that doesn’t have any ratings on Turkopticon. Actually, because of that little add-on I’ve been able to maintain a steady 98-99% approval rate ever since I began using it.

Waiting Time

So how long does it take to get judgments? I’ve restricted the available worker pool pretty strongly to ensure high quality, and it’s still only taken a few hours to get a thousand judgments.

That’s pretty awesome. I’ve worked a lot with human evaluation systems before, but always using a small in-house set of judges – and what with constraints on when those judges were available, how much they were able to work each week, and other tasks taking higher priority, it’d invariably take at least a few days before I’d receive any useful data back.

Getting thousands of judgments in a couple hours means I can launch an MTurk task when I leave for work in the morning and have it done before lunch, which makes experimenting with a lot of different ideas much faster and easier.

Scale

So how many judgments can you actually get before you run out of workers? I’m still a small fish in the MTurk system, but I’m told by my MTurk contact at Amazon that there are companies getting over a million judgments each month.

I also asked my pool of workers how much they’re available to work, in case I would need to scale up to more judgments later on, and here are some samples from what they said:

Typically, I work a total of 20-25 hours per week for a small select group of requesters. I could put in at least 20 hours per week for you alone if you were to make a custom qualification for me. If I know that I can continue to do exemplary work beyond 20 hours, I would be willing to put forth more hours of work. I want to make sure that you are getting the quality of work that you need.

On a day when I don’t have those other assignments, I’d guess I’m turking 5 to 7 hours a day (including weekends). I like to look for a large batch of HITs (preferably in the thousands) so that I can settle into a groove of being able to do them fairly quickly and once I find something like that I can happily settle in for several hours at a time.

I spend more time than the average person on mturk. I log on at about 5:30 AM and am constantly checking for work throughout the day. If the work is available, I will spend until 9PM working. Granted, I do have to take some breaks throughout the day to take care of my 3 year old, but for the most part, I am doing my best to earn while the hits are posted. If I take any time off, it is on the weekend (if I reach my earning goals for the week).

Of course, how much I can work varies. My main source of income is transcription for a market research company and mturk fills in my downtime. If I have an audio file from them, that gets my attention. If not, I’m on mturk. As a single mother working from home, I love the flexibility.

End

I’ll end with a couple other notes.

• How do other companies use human evaluation systems? Google and Bing use human judgments in their search metrics, though I think they use an in-house set of judges rather than Mechanical Turk. I’ve heard Aardvark and Quora used Mechanical Turk to seed answers when they first launched their sites. There’s also a nice set of case studies here (search for the “On-Demand Workforce” section); in particular, Knewton’s use of MTurk for performance and QA testing is pretty interesting.
• I’ve described one way of finding good workers, namely, using the filters Amazon provides. Another way could be to build a reputation system yourself, perhaps using an EM-style algorithm to determine judge quality.
• Crowdflower is another crowdsourcing system. There are a couple differences with MTurk:
  • Crowdflower’s worker pool consists of about 20 different sources, including Mechanical Turk, as well as sources like TrialPay (people can opt to complete a MTurk task to receive some kind of TrialPay deal).
  • Crowdflower offers both a self-serve platform (like MTurk), as well as a more enterprise-centric solution (where you work directly with a Crowdflower employee). The enterprise offering is pretty nice, since that means Crowdflower will take care of the lower-level details for you (like actually designing and creating the job), and they can offer suggestions for designing the HIT based on their experience.
  • Crowdflower provides the option of adding gold standard judgments to your task (items where you provide a golden answer, which are then randomly shown to workers; these are then used to monitor judges) and they try to automatically determine judge quality and item accuracy for you (e.g., by having each item judged by three different workers).
• An excellent crowdsourcing resource is CrowdScope. I also like the Deneme blog (though it hasn’t been updated in a while) for a lot of fun experiments. Panos Ipeirotis’ blog has good information as well.

Imagine you’re a budding chef. A data-curious one, of course, so you start by taking a set of foods (pizza, salad, spaghetti, etc.) and ask 10 friends how much of each they ate in the past day.

Your goal: to find natural groups of foodies, so that you can better cater to each cluster’s tastes. For example, your fratboy friends might love wings and beer, your anime friends might love soba and sushi, your hipster friends probably dig tofu, and so on.

So how can you use the data you’ve gathered to discover different kinds of groups?

One way is to use a standard clustering algorithm like k-means or Gaussian mixture modeling (see this previous post for a brief introduction). The problem is that these both assume a fixed number of clusters, which they need to be told to find. There are a couple methods for selecting the number of clusters to learn (e.g., the gap and prediction strength statistics), but the problem is a more fundamental one: most real-world data simply doesn’t have a fixed number of clusters.

That is, suppose we’ve asked 10 of our friends what they ate in the past day, and we want to find groups of eating preferences. There’s really an infinite number of foodie types (carnivore, vegan, snacker, Italian, healthy, fast food, heavy eaters, light eaters, and so on), but with only 10 friends, we simply don’t have enough data to detect them all. (Indeed, we’re limited to 10 clusters!) So whereas k-means starts with the incorrect assumption that there’s a fixed, finite number of clusters that our points come from, no matter if we feed it more data, what we’d really like is a method positing an infinite number of hidden clusters that naturally arise as we ask more friends about their food habits. (For example, with only 2 data points, we might not be able to tell the difference between vegans and vegetarians, but with 200 data points, we probably could.)

Luckily for us, this is precisely the purview of nonparametric Bayes.*

*Nonparametric Bayes refers to a class of techniques that allow some parameters to change with the data. In our case, for example, instead of fixing the number of clusters to be discovered, we allow it to grow as more data comes in.

A Generative Story

Let’s describe a generative model for finding clusters in any set of data. We assume an infinite set of latent groups, where each group is described by some set of parameters. For example, each group could be a Gaussian with a specified mean $\mu_i$ and standard deviation $\sigma_i$, and these group parameters themselves are assumed to come from some base distribution $G_0$. Data is then generated in the following manner:

• Select a cluster.
• Sample from that cluster to generate a new point.

(Note the resemblance to a finite mixture model.)

For example, suppose we ask 10 friends how many calories of pizza, salad, and rice they ate yesterday. Our groups could be:

• A Gaussian centered at (pizza = 5000, salad = 100, rice = 500) (i.e., a pizza lovers group).
• A Gaussian centered at (pizza = 100, salad = 3000, rice = 1000) (maybe a vegan group).
• A Gaussian centered at (pizza = 100, salad = 100, rice = 10000) (definitely Asian).
• …

When deciding what to eat when she woke up yesterday, Alice could have thought girl, I’m in the mood for pizza and her food consumption yesterday would have been a sample from the pizza Gaussian. Similarly, Bob could have spent the day in Chinatown, thereby sampling from the Asian Gaussian for his day’s meals. And so on.

The big question, then, is: how do we assign each friend to a group?

Assigning Groups

Chinese Restaurant Process

One way to assign friends to groups is to use a Chinese Restaurant Process. This works as follows: Imagine a restaurant where all your friends went to eat yesterday…

• Initially the restaurant is empty.
• The first person to enter (Alice) sits down at a table (selects a group). She then orders food for the table (i.e., she selects parameters for the group); everyone else who joins the table will then be limited to eating from the food she ordered.
• The second person to enter (Bob) sits down at a table. Which table does he sit at? With probability $\alpha / (1 + \alpha)$ he sits down at a new table (i.e., selects a new group) and orders food for the table; with probability $1 / (1 + \alpha)$ he sits with Alice and eats from the food she’s already ordered (i.e., he’s in the same group as Alice).
• …
• The (n+1)-st person sits down at a new table with probability $\alpha / (n + \alpha)$, and at table k with probability $n_k / (n + \alpha)$, where $n_k$ is the number of people currently sitting at table k.

Note a couple things:

• The more people (data points) there are at a table (cluster), the more likely it is that people (new data points) will join it. In other words, our groups satisfy a rich get richer property.
• There’s always a small probability that someone joins an entirely new table (i.e., a new group is formed).
• The probability of a new group depends on $\alpha$. So we can think of $\alpha$ as a dispersion parameter that affects the dispersion of our datapoints. The lower alpha is, the more tightly clustered our data points; the higher it is, the more clusters we have in any finite set of points.

(Also notice the resemblance between table selection probabilities and a Dirichlet distribution…)

Just to summarize, given n data points, the Chinese Restaurant Process specifies a distribution over partitions (table assignments) of these points. We can also generate parameters for each partition/table from a base distribution $G_0$ (for example, each table could represent a Gaussian whose mean and standard deviation are sampled from $G_0$), though to be clear, this is not part of the CRP itself.

Code

Since code makes everything better, here’s some Ruby to simulate a CRP:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

# Generate table assignments for `num_customers` customers, according to
# a Chinese Restaurant Process with dispersion parameter `alpha`.
#
# returns an array of integer table assignments
def chinese_restaurant_process(num_customers, alpha)
 return [] if num_customers <= 0

 table_assignments = [1] # first customer sits at table 1
 next_open_table = 2 # index of the next empty table

 # Now generate table assignments for the rest of the customers.
 1.upto(num_customers - 1) do |i|
   if rand < alpha.to_f / (alpha + i)
     # Customer sits at new table.
     table_assignments << next_open_table
     next_open_table += 1
   else
     # Customer sits at an existing table.
     # He chooses which table to sit at by giving equal weight to each
     # customer already sitting at a table. 
     which_table = table_assignments[rand(table_assignments.size)]
     table_assignments << which_table
   end
 end

 table_assignments
end

And here’s some sample output:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

> chinese_restaurant_process(num_customers = 10, alpha = 1)
1, 2, 3, 4, 3, 3, 2, 1, 4, 3 # table assignments from run 1
1, 1, 1, 1, 1, 1, 2, 2, 1, 3 # table assignments from run 2
1, 2, 2, 1, 3, 3, 2, 1, 3, 4 # table assignments from run 3

> chinese_restaurant_process(num_customers = 10, alpha = 3)
1, 2, 1, 1, 3, 1, 2, 3, 4, 5
1, 2, 3, 3, 4, 3, 4, 4, 5, 5
1, 1, 2, 3, 1, 4, 4, 3, 1, 1

> chinese_restaurant_process(num_customers = 10, alpha = 5)
1, 2, 1, 3, 4, 5, 6, 7, 1, 8
1, 2, 3, 3, 4, 5, 6, 5, 6, 7
1, 2, 3, 4, 5, 6, 2, 7, 2, 1

Notice that as we increase $\alpha$, so too does the number of distinct tables increase.

Polya Urn Model

Another method for assigning friends to groups is to follow the Polya Urn Model. This is basically the same model as the Chinese Restaurant Process, just with a different metaphor.

• We start with an urn containing $\alpha G_0(x)$ balls of “color” $x$, for each possible value of $x$. ($G_0$ is our base distribution, and $G_0(x)$ is the probability of sampling $x$ from $G_0$). Note that these are possibly fractional balls.
• At each time step, draw a ball from the urn, note its color, and then drop both the original ball plus a new ball of the same color back into the urn.

Note the connection between this process and the CRP: balls correspond to people (i.e., data points), colors correspond to table assignments (i.e., clusters), alpha is again a dispersion parameter (put differently, a prior), colors satisfy a rich-get-richer property (since colors with many balls are more likely to get drawn), and so on. (Again, there’s also a connection between this urn model and the urn model for the (finite) Dirichlet distribution…)

To be precise, the difference between the CRP and the Polya Urn Model is that the CRP specifies only a distribution over partitions (i.e., table assignments), but doesn’t assign parameters to each group, whereas the Polya Urn Model does both.

Code

Again, here’s some code for simulating a Polya Urn Model:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# Draw `num_balls` colored balls according to a Polya Urn Model
# with a specified base color distribution and dispersion parameter
# `alpha`.
#
# returns an array of ball colors
def polya_urn_model(base_color_distribution, num_balls, alpha)
  return [] if num_balls <= 0

  balls_in_urn = []
  0.upto(num_balls - 1) do |i|
    if rand < alpha.to_f / (alpha + balls_in_urn.size)
      # Draw a new color, put a ball of this color in the urn.
      new_color = base_color_distribution.call
      balls_in_urn << new_color
    else
      # Draw a ball from the urn, add another ball of the same color.
      ball = balls_in_urn[rand(balls_in_urn.size)]
      balls_in_urn << ball
    end
  end

  balls_in_urn
end

And here’s some sample output, using a uniform distribution over the unit interval as the color distribution to sample from:

1
2
3
4
5
6

> unit_uniform = lambda { (rand * 100).to_i / 100.0 }

> polya_urn_model(unit_uniform, num_balls = 10, alpha = 1)
0.27, 0.89, 0.89, 0.89, 0.73, 0.98, 0.43, 0.98, 0.89, 0.53 # colors in the urn from run 1
0.26, 0.26, 0.46, 0.26, 0.26, 0.26, 0.26, 0.26, 0.26, 0.85 # colors in the urn from run 2
0.96, 0.87, 0.96, 0.87, 0.96, 0.96, 0.87, 0.96, 0.96, 0.96 # colors in the urn from run 3

Code, Take 2

Here’s the same code for a Polya Urn Model, but in R:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# Return a vector of `num_balls` ball colors according to a Polya Urn Model
# with dispersion `alpha`, sampling from a specified base color distribution.
polya_urn_model = function(base_color_distribution, num_balls, alpha) {
  balls = c()

  for (i in 1:num_balls) {
    if (runif(1) < alpha / (alpha + length(balls))) {
      # Add a new ball color.
      new_color = base_color_distribution()
      balls = c(balls, new_color)
    } else {
      # Pick out a ball from the urn, and add back a
      # ball of the same color.
      ball = balls[sample(1:length(balls), 1)]
      balls = c(balls, ball)
    }
  }

  balls
}

Here are some sample density plots of the colors in the urn, when using a unit normal as the base color distribution:

Notice that as alpha increases (i.e., we sample more new ball colors from our base; i.e., as we place more weight on our prior), the colors in the urn tend to a unit normal (our base color distribution).

And here are some sample plots of points generated by the urn, for varying values of alpha:

• Each color in the urn is sampled from a uniform distribution over [0,10]x[0,10] (i.e., a [0, 10] square).
• Each group is a Gaussian with standard deviation 0.1 and mean equal to its associated color, and these Gaussian groups generate points.

Notice that the points clump together in fewer clusters for low values of alpha, but become more dispersed as alpha increases.

Stick-Breaking Process

Imagine running either the Chinese Restaurant Process or the Polya Urn Model without stop. For each group $i$, this gives a proportion $w_i$ of points that fall into group $i$.

So instead of running the CRP or Polya Urn model to figure out these proportions, can we simply generate them directly?

This is exactly what the Stick-Breaking Process does:

• Start with a stick of length one.
• Generate a random variable $\beta_1 \sim Beta(1, \alpha)$. By the definition of the Beta distribution, this will be a real number between 0 and 1, with expected value $1 / (1 + \alpha)$. Break off the stick at $\beta_1$; $w_1$ is then the length of the stick on the left.
• Now take the stick to the right, and generate $\beta_2 \sim Beta(1, \alpha)$. Break off the stick $\beta_2$ into the stick. Again, $w_2$ is the length of the stick to the left, i.e., $w_2 = (1 - \beta_1) \beta_2$.
• And so on.

Thus, the Stick-Breaking process is simply the CRP or Polya Urn Model from a different point of view. For example, assigning customers to table 1 according to the Chinese Restaurant Process is equivalent to assigning customers to table 1 with probability $w_1$.

Code

Here’s some R code for simulating a Stick-Breaking process:

1
2
3
4
5
6
7
8
9
10
11
12
13

# Return a vector of weights drawn from a stick-breaking process
# with dispersion `alpha`.
#
# Recall that the kth weight is
#   \beta_k = (1 - \beta_1) * (1 - \beta_2) * ... * (1 - \beta_{k-1}) * beta_k
# where each $\\beta\_i$ is drawn from a Beta distribution
#   \beta_i ~ Beta(1, \alpha)
stick_breaking_process = function(num_weights, alpha) {
  betas = rbeta(num_weights, 1, alpha)
  remaining_stick_lengths = c(1, cumprod(1 - betas))[1:num_weights]
  weights = remaining_stick_lengths * betas
  weights
}

And here’s some sample output:

Notice that for low values of alpha, the stick weights are concentrated on the first few weights (meaning our data points are concentrated on a few clusters), while the weights become more evenly dispersed as we increase alpha (meaning we posit more clusters in our data points).

Dirichlet Process

Suppose we run a Polya Urn Model several times, where we sample colors from a base distribution $G_0$. Each run produces a distribution of colors in the urn (say, 5% blue balls, 3% red balls, 2% pink balls, etc.), and the distribution will be different each time (for example, 5% blue balls in run 1, but 1% blue balls in run 2).

For example, let’s look again at the plots from above, where I generated samples from a Polya Urn Model with the standard unit normal as the base distribution:

Each run of the Polya Urn Model produces a slighly different distribution, though each is “centered” in some fashion around the standard Gaussian I used as base. In other words, the Polya Urn Model gives us a distribution over distributions (we get a distribution of ball colors, and this distribution of colors changes each time) – and so we finally get to the Dirichlet Process.

Formally, given a base distribution $G_0$ and a dispersion parameter $\alpha$, a sample from the Dirichlet Process $DP(G_0, \alpha)$ is a distribution $G \sim DP(G_0, \alpha)$. This sample $G$ can be thought of as a distribution of colors in a single simulation of the Polya Urn Model; sampling from $G$ gives us the balls in the urn.

So here’s the connection between the Chinese Restaurant Process, the Polya Urn Model, the Stick-Breaking Process, and the Dirichlet Process:

• Dirichlet Process: Suppose we want samples $x_i \sim G$, where $G$ is a distribution sampled from the Dirichlet Process $G \sim DP(G_0, \alpha)$.
• Polya Urn Model: One way to generate these values $x_i$ would be to take a Polya Urn Model with color distribution $G_0$ and dispersion $\alpha$. ($x_i$ would be the color of the ith ball in the urn.)
• Chinese Restaurant Process: Another way to generate $x_i$ would be to first assign tables to customers according to a Chinese Restaurant Process with dispersion $\alpha$. Every customer at the nth table would then be given the same value (color) sampled from $G_0$. ($x_i$ would be the value given to the ith customer; $x_i$ can also be thought of as the food at table $i$, or as the parameters of table $i$.)
• Stick-Breaking Process: Finally, we could generate weights $w_k$ according to a Stick-Breaking Process with dispersion $\alpha$. Next, we would give each weight $w_k$ a value (or color) $v_k$ sampled from $G_0$. Finally, we would assign $x_i$ to value (color) $v_k$ with probability $w_k$.

Recap

Let’s summarize what we’ve discussed so far.

We have a bunch of data points $p_i$ that we want to cluster, and we’ve described four essentially equivalent generative models that allow us to describe how each cluster and point could have arisen.

In the Chinese Restaurant Process:

• We generate table assignments $g_1, \ldots, g_n \sim CRP(\alpha)$ according to a Chinese Restaurant Process. ($g_i$ is the table assigned to datapoint $i$.)
• We generate table parameters $\phi_1, \ldots, \phi_n \sim G_0$ according to the base distribution $G_0$, where $\phi_k$ is the parameter for the kth distinct group.
• Given table assignments and table parameters, we generate each datapoint $p_i \sim F(\phi_{g_i})$ from a distribution $F$ with the specified table parameters. (For example, $F$ could be a Gaussian, and $\phi_i$ could be a parameter vector specifying the mean and standard deviation).

In the Polya Urn Model:

• We generate colors $\phi_1, \ldots, \phi_n \sim Polya(G_0, \alpha)$ according to a Polya Urn Model. ($\phi_i$ is the color of the ith ball.)
• Given ball colors, we generate each datapoint $p_i \sim F(\phi_i)$.

In the Stick-Breaking Process:

• We generate group probabilities (stick lengths) $w_1, \ldots, w_{\infty} \sim Stick(\alpha)$ according to a Stick-Breaking process.
• We generate group parameters $\phi_1, \ldots, \phi_{\infty} \sim G_0$ from $G_0$, where $\phi_k$ is the parameter for the kth distinct group.
• We generate group assignments $g_1, \ldots, g_n \sim Multinomial(w_1, \ldots, w_{\infty})$ for each datapoint.
• Given group assignments and group parameters, we generate each datapoint $p_i \sim F(\phi_{g_i})$.

In the Dirichlet Process:

• We generate a distribution $G \sim DP(G_0, \alpha)$ from a Dirichlet Process with base distribution $G_0$ and dispersion parameter $\alpha$.
• We generate group-level parameters $x_i \sim G$ from $G$, where $x_i$ is the group parameter for the ith datapoint. (Note: this is not the same as $\phi_i$. $x_i$ is the parameter associated to the group that the ith datapoint belongs to, whereas $\phi_k$ is the parameter of the kth distinct group.)
• Given group-level parameters $x_i$, we generate each datapoint $p_i \sim F(x_i)$.

Also, remember that each model naturally allows the number of clusters to grow as more points come in.

Inference in the Dirichlet Process Mixture

So we’ve described a generative model that allows us to calculate the probability of any particular set of group assignments to data points, but we haven’t described how to actually learn a good set of group assignments.

Let’s briefly do this now. Very roughly, the Gibbs sampling approach works as follows:

• Take the set of data points, and randomly initialize group assignments.
• Pick a point. Fix the group assignments of all the other points, and assign the chosen point a new group (which can be either an existing cluster or a new cluster) with a CRP-ish probability (as described in the models above) that depends on the group assignments and values of all the other points.
• We will eventually converge on a good set of group assignments, so repeat the previous step until happy.

For more details, this paper provides a good description. Philip Resnick and Eric Hardisty also have a friendlier, more general description of Gibbs sampling (plus an application to naive Bayes) here.

Fast Food Application: Clustering the McDonald’s Menu

Finally, let’s show an application of the Dirichlet Process Mixture. Unfortunately, I didn’t have a data set of people’s food habits offhand, so instead I took this list of McDonald’s foods and nutrition facts.

After normalizing each item to have an equal number of calories, and representing each item as a vector of (total fat, cholesterol, sodium, dietary fiber, sugars, protein, vitamin A, vitamin C, calcium, iron, calories from fat, satured fat, trans fat, carbohydrates), I ran scikit-learn’s Dirichlet Process Gaussian Mixture Model to cluster McDonald’s menu based on nutritional value.

First, how does the number of clusters inferred by the Dirichlet Process mixture vary as we feed in more (randomly ordered) points?

As expected, the Dirichlet Process model discovers more and more clusters as more and more food items arrive. (And indeed, the number of clusters appears to grow logarithmically, which can in fact be proved.)

How many clusters does the mixture model infer from the entire dataset? Running the Gibbs sampler several times, we find that the number of clusters tends around 11:

Let’s dive into one of these clusterings.

Cluster 1 (Desserts)

Looking at a sample of foods from the first cluster, we find a lot of desserts and dessert-y drinks:

• Caramel Mocha
• Frappe Caramel
• Iced Hazelnut Latte
• Iced Coffee
• Strawberry Triple Thick Shake
• Snack Size McFlurry
• Hot Caramel Sundae
• Baked Hot Apple Pie
• Cinnamon Melts
• Kiddie Cone
• Strawberry Sundae

We can also look at the nutritional profile of some foods from this cluster (after z-scaling each nutrition dimension to have mean 0 and standard deviation 1):

We see that foods in this cluster tend to be high in trans fat and low in vitamins, protein, fiber, and sodium.

Cluster 2 (Sauces)

Here’s a sample from the second cluster, which contains a lot of sauces:

• Hot Mustard Sauce
• Spicy Buffalo Sauce
• Newman’s Own Low Fat Balsamic Vinaigrette

And looking at the nutritional profile of points in this cluster, we see that it’s heavy in sodium and fat:

Cluster 3 (Burgers, Crispy Foods, High-Cholesterol)

The third cluster is very burgery:

• Hamburger
• Cheeseburger
• Filet-O-Fish
• Quarter Pounder with Cheese
• Premium Grilled Chicken Club Sandwich
• Ranch Snack Wrap
• Premium Asian Salad with Crispy Chicken
• Butter Garlic Croutons
• Sausage McMuffin
• Sausage McGriddles

It’s also high in fat and sodium, and low in carbs and sugar

Cluster 4 (Creamy Sauces)

Interestingly, even though we already found a cluster of sauces above, we discover another one as well. These sauces appear to be much more cream-based:

• Creamy Ranch Sauce
• Newman’s Own Creamy Caesar Dressing
• Coffee Cream
• Iced Coffee with Sugar Free Vanilla Syrup

Nutritionally, these sauces are higher in calories from fat, and much lower in sodium:

Cluster 5 (Salads)

Here’s a salad cluster. A lot of salads also appeared in the third cluster (along with hamburgers and McMuffins), but that’s because those salads also all contained crispy chicken. The salads in this cluster are either crisp-free or have their chicken grilled instead:

• Premium Southwest Salad with Grilled Chicken
• Premium Caesar Salad with Grilled Chicken
• Side Salad
• Premium Asian Salad without Chicken
• Premium Bacon Ranch Salad without Chicken

This is reflected in the higher content of iron, vitamin A, and fiber:

Cluster 6 (More Sauces)

Again, we find another cluster of sauces:

• Ketchup Packet
• Barbeque Sauce
• Chipotle Barbeque Sauce

These are still high in sodium, but much lower in fat compared to the other sauce clusters:

Cluster 7 (Fruit and Maple Oatmeal)

Amusingly, fruit and maple oatmeal is in a cluster by itself:

• Fruit & Maple Oatmeal

Cluster 8 (Sugary Drinks)

We also get a cluster of sugary drinks:

• Strawberry Banana Smoothie
• Wild Berry Smoothie
• Iced Nonfat Vanilla Latte
• Nonfat Hazelnut
• Nonfat Vanilla Cappuccino
• Nonfat Caramel Cappuccino
• Sweet Tea
• Frozen Strawberry Lemonade
• Coca-Cola
• Minute Maid Orange Juice

In addition to high sugar content, this cluster is also high in carbohydrates and calcium, and low in fat.

Cluster 9 (Breakfast Foods)

Here’s a cluster of high-cholesterol breakfast foods:

• Sausage McMuffin with Egg
• Sausage Burrito
• Egg McMuffin
• Bacon, Egg & Chees Biscuit
• McSkillet Burrito with Sausage
• Big Breakfast with Hotcakes

Cluster 10 (Coffee Drinks)

We find a group of coffee drinks next:

• Nonfat Cappuccino
• Nonfat Latte
• Nonfat Latte with Sugar Free Vanilla Syrup
• Iced Nonfat Latte

These are much higher in calcium and protein, and lower in sugar, than the other drink cluster above:

Cluster 11 (Apples)

Here’s a cluster of apples:

• Apple Dippers with Low Fat Caramel Dip
• Apple Slices

Vitamin C, check.

And finally, here’s an overview of all the clusters at once (using a different clustering run):

No More!

I’ll end with a couple notes:

• Kevin Knight has a hilarious introduction to Bayesian inference that describes some applications of nonparametric Bayesian techniques to computational linguistics (though I don’t think he ever quite says “nonparametric Bayes” directly).
• In the Chinese Restaurant Process, each customer sits at a single table. The Indian Buffet Process is an extension that allows customers to sample food from multiple tables (i.e., belong to multiple clusters).
• The Chinese Restaurant Process, the Polya Urn Model, and the Stick-Breaking Process are all sequential models for generating groups: to figure out table parameters in the CRP, for example, you wait for customer 1 to come in, then customer 2, then customer 3, and so on. The equivalent Dirichlet Process, on the other hand, is a parallel model for generating groups: just sample $G \sim DP(G_0, alpha)$, and then all your group parameters can be independently generated by sampling from $G$ at once. This duality is an instance of a more general phenomenon known as de Finetti’s theorem.

And that’s it.

It’s often easier to understand a chart than a table. So why is it still so hard to make a simple data graphic, and why am I still bombarded by mind-numbing reams of raw numbers?

(Yeah, I love ggplot2 to death. But sometimes I want a little more interaction, and sometimes all I want is to drag-and-drop and be done.)

So I’ve been experimenting with a small, ggplot2-inspired d3 app.

Simply drop a file, and bam! Instant scatterplot:

But wait – that’s only 2 dimensions. You can add some more through color, size, and groups:

(Click here to play with the data yourself.)

And you can easily switch which variables are getting plotted, and see all the information associated with each point.

(Same dataset, different aesthetic assignments.)

I’m thinking of adding more kinds of charts, support for categorical variables, more interactivity (sliders to interact with other dimensions?!), and making the UI even easier (e.g., simplify column naming). In the meantime, the code is here on Github, and tips and suggestions are welcome!

Scalding is an in-house MapReduce framework that Twitter recently open-sourced. Like Pig, it provides an abstraction on top of MapReduce that makes it easy to write big data jobs in a syntax that’s simple and concise. Unlike Pig, Scalding is written in pure Scala – which means all the power of Scala and the JVM is already built-in. No more UDFs, folks!

This is going to be an in-your-face introduction to Scalding, Twitter’s (Scala + Cascading) MapReduce framework.

In 140: instead of forcing you to write raw map and reduce functions, Scalding allows you to write natural code like

1
2

// Create a histogram of tweet lengths.
tweets.map('tweet -> 'length) { tweet : String => tweet.size }.groupBy('length) { _.size }

Not much different from the Ruby you’d write to compute tweet distributions over small data? Exactly.

Two notes before we begin:

• This Github repository contains all the code used.
• For a gentler introduction to Scalding, see this Getting Started guide on the Scalding wiki.

Movie Similarities

Imagine you run an online movie business, and you want to generate movie recommendations. You have a rating system (people can rate movies with 1 to 5 stars), and we’ll assume for simplicity that all of the ratings are stored in a TSV file somewhere.

Let’s start by reading the ratings into a Scalding job.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

/**
 * The input is a TSV file with three columns: (user, movie, rating).
 */
val INPUT_FILENAME = "data/ratings.tsv"

/**
 * Read in the input and give each field a type and name.
 */
val ratings = Tsv(INPUT_FILENAME, ('user, 'movie, 'rating))

/**
 * Let's also keep track of the total number of people who rated each movie.
 */
val numRaters =
  ratings
    // Put the number of people who rated each movie into a field called "numRaters".    
    .groupBy('movie) { _.size }.rename('size -> 'numRaters)

// Merge `ratings` with `numRaters`, by joining on their movie fields.
val ratingsWithSize =
  ratings.joinWithSmaller('movie -> 'movie, numRaters)

// ratingsWithSize now contains the following fields: (user, movie, rating, numRaters).

You want to calculate how similar pairs of movies are, so that if someone watches The Lion King, you can recommend films like Toy Story. So how should you define the similarity between two movies?

One way is to use their correlation:

• For every pair of movies A and B, find all the people who rated both A and B.
• Use these ratings to form a Movie A vector and a Movie B vector.
• Calculate the correlation between these two vectors.
• Whenever someone watches a movie, you can then recommend the movies most correlated with it.

Let’s start with the first two steps.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

/**
 * To get all pairs of co-rated movies, we'll join `ratings` against itself.
 * So first make a dummy copy of the ratings that we can join against.
 */
val ratings2 =
  ratingsWithSize
    .rename(('user, 'movie, 'rating, 'numRaters) -> ('user2, 'movie2, 'rating2, 'numRaters2))

/**
 * Now find all pairs of co-rated movies (pairs of movies that a user has rated) by
 * joining the duplicate rating streams on their user fields, 
 */
val ratingPairs =
  ratingsWithSize
    .joinWithSmaller('user -> 'user2, ratings2)
    // De-dupe so that we don't calculate similarity of both (A, B) and (B, A).
    .filter('movie, 'movie2) { movies : (String, String) => movies._1 < movies._2 }
    .project('movie, 'rating, 'numRaters, 'movie2, 'rating2, 'numRaters2)

// By grouping on ('movie, 'movie2), we can now get all the people who rated any pair of movies.

Before using these rating pairs to calculate correlation, let’s stop for a bit.

Since we’re explicitly thinking of movies as vectors of ratings, it’s natural to compute some very vector-y things like norms and dot products, as well as the length of each vector and the sum over all elements in each vector. So let’s compute these:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

/**
 * Compute dot products, norms, sums, and sizes of the rating vectors.
 */
val vectorCalcs =
  ratingPairs
    // Compute (x*y, x^2, y^2), which we need for dot products and norms.
    .map(('rating, 'rating2) -> ('ratingProd, 'ratingSq, 'rating2Sq)) {
      ratings : (Double, Double) =>
      (ratings._1 * ratings._2, math.pow(ratings._1, 2), math.pow(ratings._2, 2))
    }
    .groupBy('movie, 'movie2) { group =>
        group.size // length of each vector
        .sum('ratingProd -> 'dotProduct)
        .sum('rating -> 'ratingSum)
        .sum('rating2 -> 'rating2Sum)
        .sum('ratingSq -> 'ratingNormSq)
        .sum('rating2Sq -> 'rating2NormSq)
        .max('numRaters) // Just an easy way to make sure the numRaters field stays.
        .max('numRaters2)                
        // All of these operations chain together like in a builder object.
    }

To summarize, each row in vectorCalcs now contains the following fields:

• movie, movie2
• numRaters, numRaters2: the total number of people who rated each movie
• size: the number of people who rated both movie and movie2
• dotProduct: dot product between the movie vector (a vector of ratings) and the movie2 vector (also a vector of ratings)
• ratingSum, rating2sum: sum over all elements in each ratings vector
• ratingNormSq, rating2Normsq: squared norm of each vector

So let’s go back to calculating the correlation between movie and movie2. We could, of course, calculate correlation in the standard way: find the covariance between the movie and movie2 ratings, and divide by their standard deviations.

But recall that we can also write correlation in the following form:

$Corr(X, Y) = \frac{n \sum xy - \sum x \sum y}{\sqrt{n \sum x^2 - (\sum x)^2} \sqrt{n \sum y^2 - (\sum y)^2}}$

(See the Wikipedia page on correlation.)

Notice that every one of the elements in this formula is a field in vectorCalcs! So instead of using the standard calculation, we can use this form instead:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

val correlations =
  vectorCalcs
    .map(('size, 'dotProduct, 'ratingSum, 'rating2Sum, 'ratingNormSq, 'rating2NormSq) -> 'correlation) {
      val fields : (Double, Double, Double, Double, Double, Double) =>
      correlation(fields._1, fields._2, fields._3, fields._4, fields._5, fields._6)
    }

def correlation(size : Double, dotProduct : Double, ratingSum : Double,
  rating2Sum : Double, ratingNormSq : Double, rating2NormSq : Double) = {

  val numerator = size * dotProduct - ratingSum * rating2Sum
  val denominator = math.sqrt(size * ratingNormSq - ratingSum * ratingSum) * math.sqrt(size * rating2NormSq - rating2Sum * rating2Sum)

  numerator / denominator
}

And that’s it! To see the full code, check out the Github repository here.

Book Similarities

Let’s run this code over some real data. Unfortunately, I didn’t have a clean source of movie ratings available, so instead I used this dataset of 1 million book ratings.

I ran a quick command, using the handy scald.rb script that Scalding provides…

1
2

# Send the job off to a Hadoop cluster
scald.rb MovieSimilarities.scala --input ratings.tsv --output similarities.tsv

…and here’s a sample of the top output I got:

As we’d expect, we see that

• Harry Potter books are similar to other Harry Potter books
• Lord of the Rings books are similar to other Lord of the Rings books
• Tom Clancy is similar to John Grisham
• Chick lit (Summer Sisters, by Judy Blume) is similar to chick lit (Bridget Jones)

Just for fun, let’s also look at books similar to The Great Gatsby:

(Schoolboy memories, exactly.)

More Similarity Measures

Of course, there are lots of other similarity measures we could use besides correlation.

Cosine Similarity

Cosine similarity is a another common vector-based similarity measure.

1
2
3

def cosineSimilarity(dotProduct : Double, ratingNorm : Double, rating2Norm : Double) = {
  dotProduct / (ratingNorm * rating2Norm)
}

Correlation, Take II

We can also also add a regularized correlation, by (say) adding N virtual movie pairs that have zero correlation. This helps avoid noise if some movie pairs have very few raters in common (for example, The Great Gatsby had an unlikely raw correlation of 1 with many other books, due simply to the fact that those book pairs had very few ratings).

1
2
3
4
5
6
7
8
9

def regularizedCorrelation(size : Double, dotProduct : Double, ratingSum : Double,
  rating2Sum : Double, ratingNormSq : Double, rating2NormSq : Double,
  virtualCount : Double, priorCorrelation : Double) = {

  val unregularizedCorrelation = correlation(size, dotProduct, ratingSum, rating2Sum, ratingNormSq, rating2NormSq)
  val w = size / (size + virtualCount)

  w * unregularizedCorrelation + (1 - w) * priorCorrelation
}

Jaccard Similarity

Recall that one of the lessons of the Netflix prize was that implicit data can be quite useful – the mere fact that you rate a James Bond movie, even if you rate it quite horribly, suggests that you’d probably be interested in similar action films. So we can also ignore the value itself of each rating and use a set-based similarity measure like Jaccard similarity.

1
2
3
4

def jaccardSimilarity(usersInCommon : Double, totalUsers1 : Double, totalUsers2 : Double) = {
  val union = totalUsers1 + totalUsers2 - usersInCommon
  usersInCommon / union
}

Incorporation

Finally, let’s add all these similarity measures to our output.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

val PRIOR_COUNT = 10
val PRIOR_CORRELATION = 0

val similarities =
  vectorCalcs
    .map(('size, 'dotProduct, 'ratingSum, 'rating2Sum, 'ratingNormSq, 'rating2NormSq, 'numRaters, 'numRaters2) ->
      ('correlation, 'regularizedCorrelation, 'cosineSimilarity, 'jaccardSimilarity)) {

      fields : (Double, Double, Double, Double, Double, Double, Double, Double) =>

      val (size, dotProduct, ratingSum, rating2Sum, ratingNormSq, rating2NormSq, numRaters, numRaters2) = fields

      val corr = correlation(size, dotProduct, ratingSum, rating2Sum, ratingNormSq, rating2NormSq)
      val regCorr = regularizedCorrelation(size, dotProduct, ratingSum, rating2Sum, ratingNormSq, rating2NormSq, PRIOR_COUNT, PRIOR_CORRELATION)
      val cosSim = cosineSimilarity(dotProduct, math.sqrt(ratingNormSq), math.sqrt(rating2NormSq))
      val jaccard = jaccardSimilarity(size, numRaters, numRaters2)

      (corr, regCorr, cosSim, jaccard)
    }

Book Similarities Revisited

Let’s take another look at the book similarities above, now that we have these new fields.

Here are some of the top Book-Crossing pairs, sorted by their shrunk correlation:

Notice how regularization affects things: the Dark Tower pair has a pretty high raw correlation, but relatively few ratings (reducing our confidence in the raw correlation), so it ends up below the others.

And here are books similar to The Great Gatsby, this time ordered by cosine similarity:

Input Abstraction

So our code right now is tied to our specific ratings.tsv input. But what if we change the way we store our ratings, or what if we want to generate similarities for something entirely different?

Let’s abstract away our input. We’ll create a VectorSimilarities class that represents input data in the following format:

1
2
3
4
5
6
7

// This is an abstract method that returns a Pipe (aka, a stream of rating tuples).
// It takes in three symbols that name the user, item, and rating fields.
def input(userField : Symbol, itemField : Symbol, ratingField : Symbol) : Pipe

val ratings = input('user, 'item, 'rating)
// ...
// The rest of the code remains essentially the same.

Whenever we want to define a new input format, we simply subclass VectorSimilarities and provide a concrete implementation of the input method.

Book-Crossings

For example, here’s a class I could have used to generate the book recommendations above:

1
2
3
4
5
6
7
8
9
10

class BookCrossing(args : Args) extends VectorSimilarities(args) {
  override def input(userField : Symbol, itemField : Symbol, ratingField : Symbol) : Pipe = {
    val bookCrossingRatings =
      Tsv("book-crossing-ratings.tsv")
        .read
        .mapTo((0, 1, 2) -> (userField, itemField, ratingField)) { fields : (String, String, Double) => fields }

    bookCrossingRatings
  }
}

The input method simply reads from a TSV file and lets the VectorSimilarities superclass do all the work. Instant recommendations, BOOM.

Song Similarities with Twitter + iTunes

But why limit ourselves to books? We do, after all, have Twitter at our fingertips…

rated Born This Way by Lady GaGa 5 stars itun.es/iSg92N #iTunes

— gggf (@GalMusic92) February 8, 2012

Since iTunes lets you send a tweet whenever you rate a song, we can use these to generate music recommendations!

Again, we create a new class that overrides the abstract input defined in VectorSimilarities…

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

class ITunes(args : Args) extends VectorSimilarities(args) {
  // Example tweet:
  // rated New Kids On the Block: Super Hits by New Kids On the Block 5 stars http://itun.es/iSg3Fc #iTunes
  val ITUNES_REGEX = """rated (.+?) by (.+?) (\d) stars .*? #iTunes""".r

  override def input(userField : Symbol, itemField : Symbol, ratingField : Symbol) : Pipe = {
    val itunesRatings =
      // This is a Twitter-internal tweet source, but you could just as easily scrape 
      // Twitter yourself and provide your own source of tweets: https://dev.twitter.com/docs
      TweetSource()
        .mapTo('userId, 'text) { s => (s.getUserId, s.getText) }
        .filter('text) { text : String => text.contains("#iTunes") }
        .flatMap('text -> ('song, 'artist, 'rating)) {
          text : String =>
          ITUNES_REGEX.findFirstMatchIn(text).map { _.subgroups }.map { l => (l(0), l(1), l(2)) }
        }
        .rename(('userId, 'song, 'rating) -> (userField, itemField, ratingField))
        .project(userField, itemField, ratingField)

    itunesRatings
  }
}

…and snap! Here are some songs you might like if you recently listened to Beyoncé:

And some recommended songs if you like Lady Gaga:

GG Pandora.

Location Similarities with Foursquare Check-ins

But what if we don’t have explicit ratings? For example, we could be a news site that wants to generate article recommendations, and maybe we only have user visits on each story.

Or what if we want to generate restaurant or tourist recommendations, when all we know is who visits each location?

I’m at Empire State Building (350 5th Ave., btwn 33rd & 34th St., New York) 4sq.com/zZ5xGd

— Simon Ackerman (@SimonAckerman) February 8, 2012

Let’s finally make Foursquare check-ins useful. (I kid, I kid.)

Instead of using an explicit rating given to us, we can simply generate a dummy rating of 1 for each check-in. Correlation doesn’t make sense any more, but we can still pay attention to a measure like Jaccard simiilarity.

So we simply create a new class that scrapes tweets for Foursquare check-in information…

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

class Foursquare(args : Args) extends VectorSimilarities(args) {
  // Example tweet: I'm at The Ambassador (673 Geary St, btw Leavenworth & Jones, San Francisco) w/ 2 others http://4sq.com/xok3rI
  // Let's limit to New York for simplicity.
  val FOURSQUARE_REGEX = """I'm at (.+?) \(.*? New York""".r

  override def input(userField : Symbol, itemField : Symbol, ratingField : Symbol) : Pipe = {
    val foursquareCheckins =
      TweetSource()
        .mapTo('userId, 'text) { s => (s.getUserId.toLong, s.getText) }
        .flatMap('text -> ('location, 'rating)) {
          text : String =>
          FOURSQUARE_REGEX.findFirstMatchIn(text).map { _.subgroups }.map { l => (l(0), 1) }
        }
        .rename(('userId, 'location, 'rating) -> (userField, itemField, ratingField))
        .unique(userField, itemField, ratingField)

    foursquareCheckins
  }
}

…and bam! Here are locations similar to the Empire State Building:

Here are places you might want to check out, if you check-in at Bergdorf Goodman:

And here’s where to go after the Statue of Liberty:

Power of Twitter, yo.

RottenTomatoes Similarities

UPDATE: I found some movie data after all…

My review for ‘How to Train Your Dragon’ on Rotten Tomatoes: 4 1/2 stars >bit.ly/xtw3d3

— Benjamin West (@BenTheWest) February 21, 2012

So let’s use RottenTomatoes tweets to recommend movies! Here’s the code for a class that searches for RottenTomatoes tweets:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

class RottenTomatoes(args : Args) extends VectorSimilarities(args) {
  /**
   * Example tweets:
   * My review for 'Hop' on Rotten Tomatoes: 1 star > http://bit.ly/AB7Tl4
   * My review for 'The Bothersome Man (Den Brysomme mannen)' on Rotten Tomatoes: 3 stars-A muddled Playtime in Paris,... http://tmto.es/AvPoO2
   */
  val ROTTENTOMATOES_REGEX = """My review for '(.+?)' on Rotten Tomatoes: (\d) star""".r

  override val MIN_NUM_RATERS = 2
  override val MAX_NUM_RATERS = 1000
  override val MIN_INTERSECTION = 2

  override def input(userField : Symbol, itemField : Symbol, ratingField : Symbol) : Pipe = {
    val rottenTomatoesRatings =
      TweetSource()
        .mapTo('userId, 'text) { s => (s.getUserId.toLong, s.getText) }
        .flatMap('text -> ('movie, 'rating)) {
          text : String =>
          ROTTENTOMATOES_REGEX.findFirstMatchIn(text).map { _.subgroups }.map { x => (x(0), x(1).toInt) }
        }
        .rename(('userId, 'movie, 'rating) -> (userField, itemField, ratingField))
        .unique(userField, itemField, ratingField)

    rottenTomatoesRatings
  }
}

And here are the most similar movies discovered:

We see that

• Lord of the Rings, Harry Potter, and Star Wars movies are similar to other Lord of the Rings, Harry Potter, and Star Wars movies
• Big science fiction blockbusters (Avatar) are similar to big science fiction blockbusters (Inception)
• People who like one Justin Timberlake movie (Bad Teacher) also like other Justin Timberlake Movies (In Time). Similarly with Michael Fassbender (A Dangerous Method, Shame)
• Art house movies (The Tree of Life) stick together (Tinker Tailor Soldier Spy)

Let’s also look at the movies with the most negative correlation:

(The more you like loud and dirty popcorn movies (Thor) and vamp romance (Twilight), the less you like arthouse? SGTM.)

Next Steps

Hopefully I gave you a taste of the awesomeness of Scalding. To learn even more:

• Check out Scalding on Github.
• Read this Getting Started Guide on the Scalding wiki.
• Run through this code-based introduction, complete with Scalding jobs that you can run in local mode.
• Browse the API reference, which also contains many code snippets illustrating different Scalding functions (e.g., map, filter, flatMap, groupBy, count, join).
• And all the code for this post is here.

Watch out for more documentation soon, and you should most definitely follow @Scalding on Twitter for updates or to ask any questions.

Mad Props

And finally, a huge shoutout to Argyris Zymnis, Avi Bryant, and Oscar Boykin, the mastermind hackers who have spent (and continue spending!) unimaginable hours making Scalding a joy to use.

@argyris, @avibryant, @posco: Thanks for it all. #awesomejobguys #loveit

This is a bare-bones introduction to ggplot2, a visualization package in R. It assumes no knowledge of R.

For a better-looking version of this post, see this Github repository, which also contains some of the example datasets I use and a literate programming version of this tutorial.

Preview

Let’s start with a preview of what ggplot2 can do.

Given Fisher’s iris data set and one simple command…

1

qplot(Sepal.Length, Petal.Length, data = iris, color = Species)

…we can produce this plot of sepal length vs. petal length, colored by species.

Installation

You can download R here. After installation, you can launch R in interactive mode by either typing R on the command line or opening the standard GUI (which should have been included in the download).

R Basics

Vectors

Vectors are a core data structure in R, and are created with c(). Elements in a vector must be of the same type.

1
2

numbers = c(23, 13, 5, 7, 31)
names = c("edwin", "alice", "bob")

Elements are indexed starting at 1, and are accessed with [] notation.

1
2

numbers[1] # 23
names[1] # edwin

Data frames

Data frames are like matrices, but with named columns of different types (similar to database tables).

1
2
3
4
5

books = data.frame(
  title = c("harry potter", "war and peace", "lord of the rings"), # column named "title"
  author = c("rowling", "tolstoy", "tolkien"),
  num_pages = c("350", "875", "500")
)

You can access columns of a data frame with $.

1
2

books$title # c("harry potter", "war and peace", "lord of the rings")
books$author[1] # "rowling"

You can also create new columns with $.

1
2
3
4

books$num_bought_today = c(10, 5, 8)
books$num_bought_yesterday = c(18, 13, 20)
  
books$total\_num\_bought = books$num_bought_today + books$num_bought_yesterday

read.table

Suppose you want to import a TSV file into R as a data frame.

tsv file without header

For example, consider the data/students.tsv file (with columns describing each student’s age, test score, and name).

1
2
3

13   100 alice
14   95  bob
13   82  eve

We can import this file into R using read.table().

1
2
3
4
5

students = read.table("data/students.tsv",
  header = F, # file does not contain a header (`F` is short for `FALSE`), so we must manually specify column names                    
  sep = "\t", # file is tab-delimited        
  col.names = c("age", "score", "name") # column names
)

We can now access the different columns in the data frame with students$age, students$score, and students$name.

csv file with header

For an example of a file in a different format, look at the data/studentsWithHeader.tsv file.

1
2
3
4

age,score,name
13,100,alice
14,95,bob
13,82,eve

Here we have the same data, but now the file is comma-delimited and contains a header. We can import this file with

1
2
3
4

students = read.table("data/students.tsv",
  sep = ",",
  header = T  # first line contains column names, so we can immediately call `students$age`        
)

(Note: there is also a read.csv function that uses sep = "," by default.)

help

There are many more options that read.table can take. For a list of these, just type help(read.table) (or ?read.table) at the prompt to access documentation.

1
2
3

# These work for other functions as well.
help(read.table)
?read.table

ggplot2

With these R basics in place, let’s dive into the ggplot2 package.

Installation

One of R’s greatest strengths is its excellent set of packages. To install a package, you can use the install.packages() function.

1

install.packages("ggplot2")

To load a package into your current R session, use library().

1

library(ggplot2)

Scatterplots with qplot()

Let’s look at how to create a scatterplot in ggplot2. We’ll use the iris data frame that’s automatically loaded into R.

What does the data frame contain? We can use the head function to look at the first few rows.

1
2
3
4
5
6
7
8
9
10

head(iris) # by default, head displays the first 6 rows. see `?head`
head(iris, n = 10) # we can also explicitly set the number of rows to display

Sepal.Length Sepal.Width Petal.Length Petal.Width Species
         5.1         3.5          1.4         0.2  setosa
         4.9         3.0          1.4         0.2  setosa
         4.7         3.2          1.3         0.2  setosa
         4.6         3.1          1.5         0.2  setosa
         5.0         3.6          1.4         0.2  setosa
         5.4         3.9          1.7         0.4  setosa

(The data frame actually contains three types of species: setosa, versicolor, and virginica.)

Let’s plot Sepal.Length against Petal.Length using ggplot2’s qplot() function.

1
2
3
4
5

qplot(Sepal.Length, Petal.Length, data = iris)
# Plot Sepal.Length vs. Petal.Length, using data from the `iris` data frame.
# * First argument `Sepal.Length` goes on the x-axis.
# * Second argument `Petal.Length` goes on the y-axis.
# * `data = iris` means to look for this data in the `iris` data frame.    

To see where each species is located in this graph, we can color each point by adding a color = Species argument.

1

qplot(Sepal.Length, Petal.Length, data = iris, color = Species) # dude!

Similarly, we can let the size of each point denote sepal width, by adding a size = Sepal.Width argument.

1
2

qplot(Sepal.Length, Petal.Length, data = iris, color = Species, size = Petal.Width)
# We see that Iris setosa flowers have the narrowest petals.

1
2

qplot(Sepal.Length, Petal.Length, data = iris, color = Species, size = Petal.Width, alpha = I(0.7))
# By setting the alpha of each point to 0.7, we reduce the effects of overplotting.

Finally, let’s fix the axis labels and add a title to the plot.

1
2
3

qplot(Sepal.Length, Petal.Length, data = iris, color = Species,
  xlab = "Sepal Length", ylab = "Petal Length",
  main = "Sepal vs. Petal Length in Fisher's Iris data")

Other common geoms

In the scatterplot examples above, we implicitly used a point geom, the default when you supply two arguments to qplot().

1
2
3

# These two invocations are equivalent.
qplot(Sepal.Length, Petal.Length, data = iris, geom = "point")
qplot(Sepal.Length, Petal.Length, data = iris)

But we can also easily use other types of geoms to create more kinds of plots.

Barcharts: geom = “bar”

1
2
3
4
5
6
7
8
9

movies = data.frame(
  director = c("spielberg", "spielberg", "spielberg", "jackson", "jackson"),
  movie = c("jaws", "avatar", "schindler's list", "lotr", "king kong"),
  minutes = c(124, 163, 195, 600, 187)
)

# Plot the number of movies each director has.
qplot(director, data = movies, geom = "bar", ylab = "# movies")
# By default, the height of each bar is simply a count.

1
2
3

# But we can also supply a different weight.
# Here the height of each bar is the total running time of the director's movies.
qplot(director, weight = minutes, data = movies, geom = "bar", ylab = "total length (min.)")

Line charts: geom = “line”

1
2

qplot(Sepal.Length, Petal.Length, data = iris, geom = "line", color = Species)
# Using a line geom doesn't really make sense here, but hey.

1
2
3
4

# `Orange` is another built-in data frame that describes the growth of orange trees.
qplot(age, circumference, data = Orange, geom = "line",
  colour = Tree,
  main = "How does orange tree circumference vary with age?")

1
2

# We can also plot both points and lines.
qplot(age, circumference, data = Orange, geom = c("point", "line"), colour = Tree)

And that’s it with what I’ll cover.

Next Steps

I skipped over a lot of aspects of R and ggplot2 in this intro.

For example,

• There are many geoms (and other functionalities) in ggplot2 that I didn’t cover, e.g., boxplots and histograms.
• I didn’t talk about ggplot2’s layering system, or the grammar of graphics it’s based on.

So I’ll end with some additional resources on R and ggplot2.

• I don’t use it myself, but RStudio is a popular IDE for R.
• The official ggplot2 documentation is great and has lots of examples. There’s also an excellent book.
• plyr is another fantastic R package that’s also by Hadley Wickham (the author of ggplot2).
• The official R introduction is okay, but definitely not great. I haven’t found any R tutorials I really like, but I’ve heard good things about The Art of R Programming.

Introduction

Imagine you have a sequence of snapshots from a day in Justin Bieber’s life, and you want to label each image with the activity it represents (eating, sleeping, driving, etc.). How can you do this?

One way is to ignore the sequential nature of the snapshots, and build a per-image classifier. For example, given a month’s worth of labeled snapshots, you might learn that dark images taken at 6am tend to be about sleeping, images with lots of bright colors tend to be about dancing, images of cars are about driving, and so on.

By ignoring this sequential aspect, however, you lose a lot of information. For example, what happens if you see a close-up picture of a mouth – is it about singing or eating? If you know that the previous image is a picture of Justin Bieber eating or cooking, then it’s more likely this picture is about eating; if, however, the previous image contains Justin Bieber singing or dancing, then this one probably shows him singing as well.

Thus, to increase the accuracy of our labeler, we should incorporate the labels of nearby photos, and this is precisely what a conditional random field does.

Part-of-Speech Tagging

Let’s go into some more detail, using the more common example of part-of-speech tagging.

In POS tagging, the goal is to label a sentence (a sequence of words or tokens) with tags like ADJECTIVE, NOUN, PREPOSITION, VERB, ADVERB, ARTICLE.

For example, given the sentence “Bob drank coffee at Starbucks”, the labeling might be “Bob (NOUN) drank (VERB) coffee (NOUN) at (PREPOSITION) Starbucks (NOUN)”.

So let’s build a conditional random field to label sentences with their parts of speech. Just like any classifier, we’ll first need to decide on a set of feature functions $f_i$.

Feature Functions in a CRF

In a CRF, each feature function is a function that takes in as input:

• a sentence s
• the position i of a word in the sentence
• the label $l_i$ of the current word
• the label $l_{i-1}$ of the previous word

and outputs a real-valued number (though the numbers are often just either 0 or 1).

(Note: by restricting our features to depend on only the current and previous labels, rather than arbitrary labels throughout the sentence, I’m actually building the special case of a linear-chain CRF. For simplicity, I’m going to ignore general CRFs in this post.)

For example, one possible feature function could measure how much we suspect that the current word should be labeled as an adjective given that the previous word is “very”.

Features to Probabilities

Next, assign each feature function $f_j$ a weight $\lambda_j$ (I’ll talk below about how to learn these weights from the data). Given a sentence s, we can now score a labeling l of s by adding up the weighted features over all words in the sentence:

$score(l | s) = \sum_{j = 1}^m \sum_{i = 1}^n \lambda_j f_j(s, i, l_i, l_{i-1})$

(The first sum runs over each feature function $j$, and the inner sum runs over each position $i$ of the sentence.)

Finally, we can transform these scores into probabilities $p(l | s)$ between 0 and 1 by exponentiating and normalizing:

$p(l | s) = \frac{exp[score(l|s)]}{\sum_{l’} exp[score(l’|s)]} = \frac{exp[\sum_{j = 1}^m \sum_{i = 1}^n \lambda_j f_j(s, i, l_i, l_{i-1})]}{\sum_{l’} exp[\sum_{j = 1}^m \sum_{i = 1}^n \lambda_j f_j(s, i, l’_i, l’_{i-1})]}$

Example Feature Functions

So what do these feature functions look like? Examples of POS tagging features could include:

• $f_1(s, i, l_i, l_{i-1}) = 1$ if $l_i =$ ADVERB and the ith word ends in “-ly”; 0 otherwise.
  • If the weight $\lambda_1$ associated with this feature is large and positive, then this feature is essentially saying that we prefer labelings where words ending in -ly get labeled as ADVERB.
• $f_2(s, i, l_i, l_{i-1}) = 1$ if $i = 1$, $l_i =$ VERB, and the sentence ends in a question mark; 0 otherwise.
  • Again, if the weight $\lambda_2$ associated with this feature is large and positive, then labelings that assign VERB to the first word in a question (e.g., “Is this a sentence beginning with a verb?”) are preferred.
• $f_3(s, i, l_i, l_{i-1}) = 1$ if $l_{i-1} =$ ADJECTIVE and $l_i =$ NOUN; 0 otherwise.
  • Again, a positive weight for this feature means that adjectives tend to be followed by nouns.
• $f_4(s, i, l_i, l_{i-1}) = 1$ if $l_{i-1} =$ PREPOSITION and $l_i =$ PREPOSITION.
  • A negative weight $\lambda_4$ for this function would mean that prepositions don’t tend to follow prepositions, so we should avoid labelings where this happens.

And that’s it! To sum up: to build a conditional random field, you just define a bunch of feature functions (which can depend on the entire sentence, a current position, and nearby labels), assign them weights, and add them all together, transforming at the end to a probability if necessary.

Now let’s step back and compare CRFs to some other common machine learning techniques.

Smells like Logistic Regression…

The form of the CRF probabilities $p(l | s) = \frac{exp[\sum_{j = 1}^m \sum_{i = 1}^n f_j(s, i, l_i, l_{i-1})]}{\sum_{l’} exp[\sum_{j = 1}^m \sum_{i = 1}^n f_j(s, i, l’_i, l’_{i-1})]}$ might look familiar.

That’s because CRFs are indeed basically the sequential version of logistic regression: whereas logistic regression is a log-linear model for classification, CRFs are a log-linear model for sequential labels.

Looks like HMMs…

Recall that Hidden Markov Models are another model for part-of-speech tagging (and sequential labeling in general). Whereas CRFs throw any bunch of functions together to get a label score, HMMs take a generative approach to labeling, defining

$p(l,s) = p(l_1) \prod_i p(l_i | l_{i-1}) p(w_i | l_i)$

where

• $p(l_i | l_{i-1})$ are transition probabilities (e.g., the probability that a preposition is followed by a noun);
• $p(w_i | l_i)$ are emission probabilities (e.g., the probability that a noun emits the word “dad”).

So how do HMMs compare to CRFs? CRFs are more powerful – they can model everything HMMs can and more. One way of seeing this is as follows.

Note that the log of the HMM probability is $\log p(l,s) = \log p(l_0) + \sum_i \log p(l_i | l_{i-1}) + \sum_i \log p(w_i | l_i)$. This has exactly the log-linear form of a CRF if we consider these log-probabilities to be the weights associated to binary transition and emission indicator features.

That is, we can build a CRF equivalent to any HMM by…

• For each HMM transition probability $ p(l_i = y | l_{i-1} = x) $, define a set of CRF transition features of the form $f_{x,y}(s, i, l_i, l_{i-1}) = 1$ if $l_i = y$ and $l_{i-1} = x$. Give each feature a weight of $w_{x,y} = \log p(l_i = y | l_{i-1} = x)$.
• Similarly, for each HMM emission probability $p(w_i = z | l_{i} = x)$, define a set of CRF emission features of the form $g_{x,y}(s, i, l_i, l_{i-1}) = 1$ if $w_i = z$ and $l_i = x$. Give each feature a weight of $w_{x,z} = \log p(w_i = z | l_i = x)$.

Thus, the score $p(l|s)$ computed by a CRF using these feature functions is precisely proportional to the score computed by the associated HMM, and so every HMM is equivalent to some CRF.

However, CRFs can model a much richer set of label distributions as well, for two main reasons:

• CRFs can define a much larger set of features. Whereas HMMs are necessarily local in nature (because they’re constrained to binary transition and emission feature functions, which force each word to depend only on the current label and each label to depend only on the previous label), CRFs can use more global features. For example, one of the features in our POS tagger above increased the probability of labelings that tagged the first word of a sentence as a VERB if the end of the sentence contained a question mark.
• CRFs can have arbitrary weights. Whereas the probabilities of an HMM must satisfy certain constraints (e.g., $0 <= p(w_i | l_i) <= 1, \sum_w p(w_i = w | l_1) = 1)$, the weights of a CRF are unrestricted (e.g., $\log p(w_i | l_i)$ can be anything it wants).

Learning Weights

Let’s go back to the question of how to learn the feature weights in a CRF. One way is (surprise) to use gradient ascent.

Assume we have a bunch of training examples (sentences and associated part-of-speech labels). Randomly initialize the weights of our CRF model. To shift these randomly initialized weights to the correct ones, for each training example…

• Go through each feature function $f_i$, and calculate the gradient of the log probability of the training example with respect to $\lambda_i$: $\frac{\partial}{\partial w_j} \log p(l | s) = \sum_{j = 1}^m f_i(s, j, l_j, l_{j-1}) - \sum_{l’} p(l’ | s) \sum_{j = 1}^m f_i(s, j, l’_j, l’_{j-1})$
• Note that the first term in the gradient is the contribution of feature $f_i$ under the true label, and the second term in the gradient is the expected contribution of feature $f_i$ under the current model. This is exactly the form you’d expect gradient ascent to take.
• Move $\lambda_i$ in the direction of the gradient: $\lambda_i = \lambda_i + \alpha [\sum_{j = 1}^m f_i(s, j, l_j, l_{j-1}) - \sum_{l’} p(l’ | s) \sum_{j = 1}^m f_i(s, j, l’_j, l’_{j-1})]$ where $\alpha$ is some learning rate.
• Repeat the previous steps until some stopping condition is reached (e.g., the updates fall below some threshold).

In other words, every step takes the difference between what we want the model to learn and the model’s current state, and moves $\lambda_i$ in the direction of this difference.

Finding the Optimal Labeling

Suppose we’ve trained our CRF model, and now a new sentence comes in. How do we do label it?

The naive way is to calculate $p(l | s)$ for every possible labeling l, and then choose the label that maximizes this probability. However, since there are $k^m$ possible labels for a tag set of size k and a sentence of length m, this approach would have to check an exponential number of labels.

A better way is to realize that (linear-chain) CRFs satisfy an optimal substructure property that allows us to use a (polynomial-time) dynamic programming algorithm to find the optimal label, similar to the Viterbi algorithm for HMMs.

A More Interesting Application

Okay, so part-of-speech tagging is kind of boring, and there are plenty of existing POS taggers out there. When might you use a CRF in real life?

Suppose you want to mine Twitter for the types of presents people received for Christmas:

What people on Twitter wanted for Christmas, and what they got: twitter.com/edchedch/statu…

— Edwin Chen (@edchedch) January 2, 2012

(Yes, I just embedded a tweet. BOOM.)

How can you figure out which words refer to gifts?

To gather data for the graphs above, I simply looked for phrases of the form “I want XXX for Christmas” and “I got XXX for Christmas”. However, a more sophisticated CRF variant could use a GIFT part-of-speech-like tag (even adding other tags like GIFT-GIVER and GIFT-RECEIVER, to get even more information on who got what from whom) and treat this like a POS tagging problem. Features could be based around things like “this word is a GIFT if the previous word was a GIFT-RECEIVER and the word before that was ‘gave’” or “this word is a GIFT if the next two words are ‘for Christmas’”.

Fin

I’ll end with some more random thoughts:

• I explicitly skipped over the graphical models framework that conditional random fields sit in, because I don’t think they add much to an initial understanding of CRFs. But if you’re interested in learning more, Daphne Koller is teaching a free, online course on graphical models starting in January.
• Or, if you’re more interested in the many NLP applications of CRFs (like part-of-speech tagging or named entity extraction), Manning and Jurafsky are teaching an NLP class in the same spirit.
• I also glossed a bit over the analogy between CRFs:HMMs and Logistic Regression:Naive Bayes. This image (from Sutton and McCallum’s introduction to conditional random fields) sums it up, and shows the graphical model nature of CRFs as well:

How was the Netflix Prize won? I went through a lot of the Netflix Prize papers a couple years ago, so I’ll try to give an overview of the techniques that went into the winning solution here.

Normalization of Global Effects

Suppose Alice rates Inception 4 stars. We can think of this rating as composed of several parts:

• A baseline rating (e.g., maybe the mean over all user-movie ratings is 3.1 stars).
• An Alice-specific effect (e.g., maybe Alice tends to rate movies lower than the average user, so her ratings are -0.5 stars lower than we normally expect).
• An Inception-specific effect (e.g., Inception is a pretty awesome movie, so its ratings are 0.7 stars higher than we normally expect).
• A less predictable effect based on the specific interaction between Alice and Inception that accounts for the remainder of the stars (e.g., Alice really liked Inception because of its particular combination of Leonardo DiCaprio and neuroscience, so this rating gets an additional 0.7 stars).

In other words, we’ve decomposed the 4-star rating into: 4 = [3.1 (the baseline rating) - 0.5 (the Alice effect) + 0.7 (the Inception effect)] + 0.7 (the specific interaction)

So instead of having our models predict the 4-star rating itself, we could first try to remove the effect of the baseline predictors (the first three components) and have them predict the specific 0.7 stars. (I guess you can also think of this as a simple kind of boosting.)

More generally, additional baseline predictors include:

• A factor that allows Alice’s rating to (linearly) depend on the (square root of the) number of days since her first rating. (For example, have you ever noticed that you become a harsher critic over time?)
• A factor that allows Alice’s rating to depend on the number of days since the movie’s first rating by anyone. (If you’re one of the first people to watch it, maybe it’s because you’re a huge fan and really excited to see it on DVD, so you’ll tend to rate it higher.)
• A factor that allows Alice’s rating to depend on the number of people who have rated Inception. (Maybe Alice is a hipster who hates being part of the crowd.)
• A factor that allows Alice’s rating to depend on the movie’s overall rating.
• (Plus a bunch of others.)

And, in fact, modeling these biases turned out to be fairly important: in their paper describing their final solution to the Netflix Prize, Bell and Koren write that

Of the numerous new algorithmic contributions, I would like to highlight one – those humble baseline predictors (or biases), which capture main effects in the data. While the literature mostly concentrates on the more sophisticated algorithmic aspects, we have learned that an accurate treatment of main effects is probably at least as signficant as coming up with modeling breakthroughs.

(For a perhaps more concrete example of why removing these biases is useful, suppose you know that Bob likes the same kinds of movies that Alice does. To predict Bob’s rating of Inception, instead of simply predicting the same 4 stars that Alice rated, if we know that Bob tends to rate movies 0.3 stars higher than average, then we could first remove Alice’s bias and then add in Bob’s: 4 + 0.5 + 0.3 = 4.8.)

Neighborhood Models

Let’s now look at some slightly more sophisticated models. As alluded to in the section above, one of the standard approaches to collaborative filtering is to use neighborhood models.

Briefly, a neighborhood model works as follows. To predict Alice’s rating of Titanic, you could do two things:

• Item-item approach: find a set of items similar to Titanic that Alice has also rated, and take the (weighted) mean of Alice’s ratings on them.
• User-user approach: find a set of users similar to Alice who rated Titanic, and again take the mean of their ratings of Titanic.

(See also my post on item-to-item collaborative filtering on Amazon.)

The main questions, then, are (let’s stick to the item-item approach for simplicity):

• How do we find the set of similar items?
• How do we weight these items when taking their mean?

The standard approach is to take some similarity metric (e.g., correlation or a Jaccard index) to define similarities between pairs of movies, take the K most similar movies under this metric (where K is perhaps chosen via cross-validation), and then use the same similarity metric when computing the weighted mean.

This has a couple problems:

• Neighbors aren’t independent, so using a standard similarity metric to define a weighted mean overcounts information. For example, suppose you ask five friends where you should eat tonight. Three of them went to Mexico last week and are sick of burritos, so they strongly recommend against a taqueria. Thus, your friends’ recommendations have a stronger bias than what you’d get if you asked five friends who didn’t know each other at all. (Compare with the situation where all three Lord of the Rings Movies are neighbors of Harry Potter.)
• Different movies should perhaps be using different numbers of neighbors. Some movies may be predicted well by only one neighbor (e.g., Harry Potter 2 could be predicted well by Harry Potter 1 alone), some movies may require more, and some movies may have no good neighbors (so you should ignore your neighborhood algorithms entirely and let your other ratings models stand on their own).

So another approach is the following:

• You can still use a similarity metric like correlation or cosine similarity to choose the set of similar items.
• But instead of using the similarity metric to define the interpolation weights in the mean calculations, you essentially perform a (sparse) linear regression to find the weights that minimize the squared error between an item’s rating and a linear combination of the ratings of its neighbors. Note that these weights are no longer constrained, so that if all neighbors are weak, then their weights will be close to zero and the neighborhood model will have a low effect.

(A slightly more complicated user-user approach, similar to this item-item neighborhood approach, is also useful.)

Implicit Data

Adding on to the neighborhood approach, we can also let implicit data influence our predictions. The mere fact that a user rated lots of science fiction movies but no westerns, suggests that the user likes science fiction better than cowboys. So using a similar framework as in the neighborhood ratings model, we can learn for Inception a set of offset weights associated to Inception’s movie neighbors.

Whenever we want to predict how Bob rates Inception, we look at whether Bob rated each of Inception’s neighbors. If he did, we add in the corresponding offset; if not, then we add nothing (and, thus, Bob’s rating is implicitly penalized by the missing weight).

Matrix Factorization

Complementing the neighborhood approach to collaborative filtering is the matrix factorization approach. Whereas the neighborhood approach takes a very local approach to ratings (if you liked Harry Potter 1, then you’ll like Harry Potter 2!), the factorization approach takes a more global view (we know that you like fantasy movies and that Harry Potter has a strong fantasy element, so we think that you’ll like Harry Potter) that decomposes users and movies into a set of latent factors (which we can think of as categories like “fantasy” or “violence”).

In fact, matrix factorization methods were probably the most important class of techniques for winning the Netflix Prize. In their 2008 Progress Prize paper, Bell and Koren write

It seems that models based on matrix-factorization were found to be most accurate (and thus popular), as evident by recent publications and discussions on the Netflix Prize forum. We definitely agree to that, and would like to add that those matrix-factorization models also offer the important flexibility needed for modeling temporal effects and the binary view. Nonetheless, neighborhood models, which have been dominating most of the collaborative filtering literature, are still expected to be popular due to their practical characteristics - being able to handle new users/ratings without re-training and offering direct explanations to the recommendations.

The typical way to perform matrix factorizations is to perform a singular value decomposition on the (sparse) ratings matrix (using stochastic gradient descent and regularizing the weights of the factors, possibly constraining the weights to be positive to get a type of non-negative matrix factorization). (Note that this “SVD” is a little different from the standard SVD learned in linear algebra, since not every user has rated every movie and so the ratings matrix contains many missing elements that we don’t want to simply treat as 0.)

Some SVD-inspired methods used in the Netflix Prize include:

• Standard SVD: Once you’ve represented users and movies as factor vectors, you can dot product Alice’s vector with Inception’s vector to get Alice’s predicted rating of Inception.
• Asymmetric SVD: Instead of users having their own notion of factor vectors, we can represent users as a bag of items they have rated (or provided implicit feedback for). So Alice is now represented as a (possibly weighted) sum of the factor vectors of the items she has rated, and to get her predicted rating of Titanic, we can dot product this representation with the factor vector of Titanic. From a practical perspective, this model has an added benefit in that no user parameterizations are needed, so we can use this approach to generate recommendations as soon as a user provides some feedback (which could just be views or clicks on an item, and not necessarily ratings), without needing to retrain the model to factorize the user.
• SVD++: Incorporate both the standard SVD and the asymmetric SVD model by representing users both by their own factor representation and as a bag of item vectors.

Regression

Some regression models were also used in the predictions. The models are fairly standard, I think, so I won’t spend too long here. Basically, just as with the neighborhood models, we can take a user-centric approach and a movie-centric approach to regression:

• User-centric approach: We learn a regression model for each user, using all the movies that the user rated as the dataset. The response is the movie’s rating, and the predictor variables are attributes associated to that movie (which can be derived from, say, PCA, MDS, or an SVD).
• Movie-centric approach: Similarly, we can learn a regression model for each movie, using all the users that rated the movie as the dataset.

Restricted Boltzmann Machines

Restricted Boltzmann Machines provide another kind of latent factor approach that can be used. See this paper for a description of how to apply them to the Netflix Prize. (In case the paper’s a little difficult to read, I wrote an introduction to RBMs a little while ago.)

Temporal Effects

Many of the models incorporate temporal effects. For example, when describing the baseline predictors above, we used a few temporal predictors that allowed a user’s rating to (linearly) depend on the time since the first rating he ever made and on the time since a movie’s first rating. We can also get more fine-grained temporal effects by, say, binning items into a couple months’ worth of ratings at a time, and allowing movie biases to change within each bin. (For example, maybe in May 2006, Time Magazine nominated Titanic as the best movie ever made, which caused a spurt in glowing ratings around that time.)

In the matrix factorization approach, user factors were also allowed to be time-dependent (e.g., maybe Bob comes to like comedy movies more and more over time). We can also give more weight to recent user actions.

Regularization

Regularization was also applied throughout pretty much all the models learned, to prevent overfitting on the dataset. Ridge regression was heavily used in the factorization models to penalize large weights, and lasso regression (though less effective) was useful as well. Many other parameters (e.g., the baseline predictors, similarity weights and interpolation weights in the neighborhood models) were also estimated using fairly standard shrinkage techniques.

Ensemble Methods

Finally, let’s talk about how all of these different algorithms were combined to provide a single rating that exploits the strengths of each model. (Note that, as mentioned above, many of these models were not trained on the raw ratings data directly, but rather on the residuals of other models.)

In the paper detailing their final solution, the winners describe using gradient boosted decision trees to combine over 500 models; previous solutions used instead a linear regression to combine the predictors.

Briefly, gradient boosted decision trees work by sequentially fitting a series of decision trees to the data; each tree is asked to predict the error made by the previous trees, and is often trained on slightly perturbed versions of the data. (For a longer description of a similar technique, see my introduction to random forests.)

Since GBDTs have a built-in ability to apply different methods to different slices of the data, we can add in some predictors that help the trees make useful clusterings:

• Number of movies each user rated
• Number of users that rated each movie
• Factor vectors of users and movies
• Hidden units of a restricted Boltzmann Machine

(For example, one thing that Bell and Koren found (when using an earlier ensemble method) was that RBMs are more useful when the movie or the user has a low number of ratings, and that matrix factorization methods are more useful when the movie or user has a high number of ratings.)

Here’s a graph of the effect of ensemble size from early on in the competition (in 2007), and the authors’ take on it:

However, we would like to stress that it is not necessary to have such a large number of models to do well. The plot below shows RMSE as a function of the number of methods used. One can achieve our winning score (RMSE=0.8712) with less than 50 methods, using the best 3 methods can yield RMSE < 0.8800, which would land in the top 10. Even just using our single best method puts us on the leaderboard with an RMSE of 0.8890. The lesson here is that having lots of models is useful for the incremental results needed to win competitions, but practically, excellent systems can be built with just a few well-selected models.

What types of students go to which schools? There are, of course, the classic stereotypes:

• MIT has the hacker engineers.
• Stanford has the laid-back, social folks.
• Harvard has the prestigious leaders of the world.
• Berkeley has the activist hippies.
• Caltech has the hardcore science nerds.

But how well do these perceptions match reality? What are students at Stanford, Harvard, MIT, Caltech, and Berkeley really interested in? Following the path of my previous data-driven post on differences between Silicon Valley and NYC, I scraped the Quora profiles of a couple hundred followers of each school to find out.

Topics

So let’s look at what kinds of topics followers of each school are interested in*. (Skip past the lists for a discussion.)

MIT

Topics are followed by p(school = MIT|topic).

• MIT Media Lab 0.893
• Ksplice 0.69
• Lisp (programming language) 0.677
• Nokia 0.659
• Public Speaking 0.65
• Data Storage 0.65
• Google Voice 0.609
• Hacking 0.602
• Startups in Europe 0.597
• Startup Names 0.572
• Mechanical Engineering 0.563
• Engineering 0.563
• Distributed Databases 0.544
• StackOverflow 0.536
• Boston 0.513
• Learning 0.507
• Open Source 0.498
• Cambridge 0.496
• Public Relations 0.493
• Visualization 0.492
• Semantic Web 0.486
• Andreessen-Horowitz 0.483
• Nature 0.475
• Cryptography 0.474
• Startups in Boston 0.452
• Adobe Photoshop 0.451
• Computer Security 0.447
• Sachin Tendulkar 0.443
• Hacker News 0.442
• Games 0.429
• Android Applications 0.428
• Best Engineers and Programmers 0.427
• College Admissions & Getting Into College 0.422
• Co-Founders 0.419
• Big Data 0.41
• System Administration 0.4
• Biotechnology 0.398
• Higher Education 0.394
• NoSQL 0.387
• User Experience 0.386
• Career Advice 0.377
• Artificial Intelligence 0.375
• Scalability 0.37
• Taylor Swift 0.368
• Google Search 0.368
• Functional Programming 0.365
• Bing 0.363
• Bioinformatics 0.361
• How I Met Your Mother (TV series) 0.361
• Operating Systems 0.356
• Compilers 0.355
• Google Chrome 0.354
• Management & Organizational Leadership 0.35
• Literary Fiction 0.35
• Intelligence 0.348
• Fight Club (1999 movie) 0.344
• Hip Hop Music 0.34
• UX Design 0.337
• Web Application Frameworks 0.336
• Startups in New York City 0.333
• Book Recommendations 0.33
• Engineering Recruiting 0.33
• Search Engines 0.329
• Social Search 0.329
• Data Science 0.328
• History 0.328
• Interaction Design 0.326
• Classification (machine learning) 0.322
• Startup Incubators and Seed Programs 0.321
• Graphic Design 0.321
• Product Design (software) 0.319
• The College Experience 0.319
• Writing 0.319
• MapReduce 0.318
• Database Systems 0.315
• User Interfaces 0.314
• Literature 0.314
• C (programming language) 0.314
• Television 0.314
• Reading 0.313
• Usability 0.312
• Books 0.312
• Computers 0.311
• Stealth Startups 0.311
• Daft Punk 0.31
• Healthy Eating 0.309
• Innovation 0.309
• Skiing 0.305
• JavaScript 0.304
• Rock Music 0.304
• Mozilla Firefox 0.304
• Self-Improvement 0.303
• McKinsey & Company 0.302
• AngelList 0.301
• Data Visualization 0.301
• Cassandra (database) 0.301

Stanford

Topics are followed by p(school = Stanford|topic).

• Stanford Computer Science 0.951
• Stanford Graduate School of Business 0.939
• Stanford 0.896
• Stanford Football 0.896
• Stanford Cardinal 0.896
• Social Dance 0.847
• Stanford University Courses 0.847
• Romance 0.769
• Instagram 0.745
• College Football 0.665
• Mobile Location Applications 0.634
• Online Communities 0.621
• Interpersonal Relationships 0.585
• Food & Restaurants in Palo Alto 0.572
• Your 20s 0.566
• Men’s Fashion 0.548
• Flipboard 0.537
• Inception (2010 movie) 0.535
• Tumblr 0.531
• People Skills 0.522
• Exercise 0.52
• Joel Spolsky 0.516
• Valuations 0.515
• The Social Network (2010 movie) 0.513
• LeBron James 0.506
• Northern California 0.506
• Evernote 0.5
• Quora Community 0.5
• Blogging 0.49
• Downtown Palo Alto 0.487
• The College Experience 0.485
• Consumer Internet 0.477
• Restaurants in San Francisco 0.477
• Chad Hurley 0.47
• Meditation 0.468
• Yishan Wong 0.466
• Arrested Development (TV series) 0.463
• fbFund 0.457
• Best Engineers at X Company 0.451
• Language 0.45
• Words 0.448
• Happiness 0.447
• Path (company) 0.446
• Color Labs (startup) 0.446
• Palo Alto 0.445
• Woot.com 0.442
• Beer 0.442
• PayPal 0.441
• Women in Startups 0.438
• Techmeme 0.433
• Women in Engineering 0.428
• The Mission (San Francisco neighborhood) 0.427
• iPhone Applications 0.416
• Asana 0.413
• Monetization 0.412
• Repetitive Strain Injury (RSI) 0.4
• IDEO 0.398
• Spotify 0.397
• San Francisco Giants 0.396
• Fortune Magazine 0.389
• Love 0.387
• Human-Computer Interaction 0.382
• Hip Hop Music 0.378
• Self-Improvement 0.378
• Food in San Francisco 0.375
• Quora (company) 0.374
• Quora Infrastructure 0.373
• iPhone 0.371
• Square (company) 0.369
• Social Psychology 0.369
• Network Effects 0.366
• Chris Sacca 0.365
• Walt Mossberg 0.364
• Salesforce.com 0.362
• Sex 0.361
• Etiquette 0.361
• David Pogue 0.361
• Gowalla 0.36
• iOS Development 0.354
• Palantir Technologies 0.353
• Mobile Computing 0.347
• Sports 0.346
• Video Games 0.345
• Burning Man 0.345
• Engineering Management 0.343
• Cognitive Science 0.342
• Dating & Relationships 0.341
• Fred Wilson (venture investor) 0.337
• Taiwan 0.333
• Natural Language Processing 0.33
• Eric Schmidt 0.329
• Social Advice 0.329
• Engineering Recruiting 0.328
• Job Interviews 0.325
• Mobile Phones 0.324
• Twitter Inc. (company) 0.321
• Engineering in Silicon Valley 0.321
• San Francisco Bay Area 0.321
• Google Analytics 0.32
• Fashion 0.315
• Interaction Design 0.314
• Open Graph 0.313
• Drugs & Pharmaceuticals 0.312
• Electronic Music 0.312
• Facebook Inc. (company) 0.309
• Fitness 0.309
• YouTube 0.308
• TED Talks 0.308
• Freakonomics (2005 Book) 0.307
• Jack Dorsey 0.306
• Nutrition 0.305
• Puzzles 0.305
• Silicon Valley Mergers & Acquisitions 0.304
• Viral Growth & Analytics 0.304
• Amazon Web Services 0.304
• StumbleUpon 0.303
• Exceptional Comment Threads 0.303

Harvard

• Harvard Business School 0.968
• Harvard Business Review 0.922
• Harvard Square 0.912
• Harvard Law School 0.912
• Jimmy Fallon 0.899
• Boston Red Sox 0.658
• Klout 0.644
• Oprah Winfrey 0.596
• Ivanka Trump 0.587
• Dalai Lama 0.569
• Food in New York City 0.565
• U2 0.562
• TwitPic 0.534
• 37signals 0.522
• David Lynch (director) 0.512
• Al Gore 0.508
• TechStars 0.49
• Baseball 0.487
• Private Equity 0.471
• Classical Music 0.46
• Startups in New York City 0.458
• HootSuite 0.449
• Kiva 0.442
• Ultimate Frisbee 0.441
• Huffington Post 0.436
• New York City 0.433
• Charlie Cheever 0.433
• The New York Times 0.431
• Technology Journalism 0.431
• McKinsey & Company 0.427
• TweetDeck 0.422
• How Does X Work? 0.417
• Ashton Kutcher 0.414
• Coldplay 0.402
• Conan O’Brien 0.397
• Fast Company 0.397
• WikiLeaks 0.394
• Michael Jackson 0.389
• Guy Kawasaki 0.389
• Journalism 0.384
• Wall Street Journal 0.384
• Cambridge 0.371
• Seattle 0.37
• Cities & Metro Areas 0.357
• Boston 0.353
• Tim Ferriss (author) 0.35
• The New Yorker 0.343
• Law 0.34
• Mashable 0.338
• Politics 0.335
• The Economist 0.334
• Barack Obama 0.333
• Skiing 0.329
• McKinsey Quarterly 0.325
• Wired (magazine) 0.316
• Bill Gates 0.31
• Mad Men (TV series) 0.308
• India 0.306
• TED Talks 0.306
• Netflix 0.304
• Wine 0.303
• Angel Investors 0.302
• Facebook Ads 0.301

UC Berkeley

• Berkeley 0.978
• UC Riverside 0.92
• California Golden Bears 0.91
• Internships 0.717
• Web Marketing 0.484
• Google Social Strategy 0.453
• Southwest Airlines 0.451
• WordPress 0.429
• Stock Market 0.429
• BMW (automobile) 0.428
• Web Applications 0.423
• Flickr 0.422
• Snowboarding 0.42
• Electronic Music 0.404
• MySQL 0.401
• Internet Advertising 0.399
• Search Engine Optimization (SEO) 0.398
• Yelp 0.396
• Groupon 0.393
• In-N-Out Burger 0.391
• The Matrix (1999 movie) 0.389
• Trading (finance) 0.385
• jQuery 0.381
• Hedge Funds 0.378
• Social Media Marketing 0.377
• San Francisco 0.376
• Stealth Startups 0.362
• Yahoo! 0.36
• Cascading Style Sheets 0.359
• Angel Investors 0.355
• UX Design 0.35
• StarCraft 0.348
• Los Angeles Lakers 0.347
• Mountain View 0.345
• How I Met Your Mother (TV series) 0.338
• Google+ 0.337
• Ruby on Rails 0.333
• Reading 0.333
• Social Media 0.326
• China 0.322
• Palantir Technologies 0.319
• Facebook Platform 0.315
• Basketball 0.315
• Education 0.314
• Business Development 0.312
• Online & Mobile Payments 0.305
• Restaurants in San Francisco 0.302
• Technology Companies 0.302
• Seth Godin 0.3

Caltech

• Pasadena 0.969
• Chess 0.748
• Table Tennis 0.671
• UCLA 0.67
• MacBook Pro 0.618
• Physics 0.618
• Haskell 0.582
• Los Angeles 0.58
• Electrical Engineering 0.567
• Star Trek (movie 0.561
• Disruptive Technology 0.545
• Science 0.53
• Biology 0.526
• Quantum Mechanics 0.521
• LaTeX 0.514
• Mathematics 0.488
• xkcd 0.488
• Genetics & Heredity 0.487
• Chemistry 0.47
• Medicine & Healthcare 0.448
• Poker 0.445
• C++ (programming language) 0.442
• Data Structures 0.434
• Emacs 0.428
• MongoDB 0.423
• Neuroscience 0.404
• Science Fiction 0.4
• Mac OS X 0.394
• Board Games 0.387
• Computers 0.386
• Research 0.385
• Finance 0.385
• The Future 0.379
• Linux 0.378
• The Colbert Report 0.376
• The Beatles 0.374
• The Onion 0.365
• Ruby 0.363
• Cars & Automobiles 0.361
• Quantitative Finance 0.359
• Academia 0.359
• Law 0.355
• Cooking 0.354
• Psychology 0.349
• Eminem 0.347
• Football (Soccer) 0.346
• Computer Programming 0.343
• Algorithms 0.343
• Evolutionary Biology 0.337
• Behavioral Economics 0.335
• California 0.329
• Machine Learning 0.326
• Futurama 0.324
• Social Advice 0.324
• StarCraft II 0.319
• Job Interview Questions 0.318
• Game Theory 0.316
• This American Life 0.315
• Economics 0.314
• Vim 0.31
• Graduate School 0.309
• Git (revision control) 0.306
• Computer Science 0.303

What do we see?

• First, in a nice validation of this approach, we find that each school is interested in exactly the locations we’d expect: Caltech is interested in Pasadena and Los Angeles; MIT and Harvard are both interested in Boston and Cambridge (Harvard is interested in New York City as well); Stanford is interested in Palo Alto, Northern California, and San Francisco Bay Area; and Berkeley is interested in Berkeley, San Francisco, and Mountain View.
• More interestingly, let’s look at where each school likes to eat. Stereotypically, we expect Harvard, Stanford, and Berkeley students to be more outgoing and social, and MIT and Caltech students to be more introverted. This is indeed what we find:
  • Harvard follows Food in New York City; Stanford follows Food & Restaurants in Palo Alto, Restaurants in San Francisco, and Food in San Francisco; and Berkeley follows Restaurants in San Francisco and In-N-Out Burger. In other words, Harvard, Stanford, and Berkeley love eating out.
  • Caltech, on the other hand, loves Cooking, and MIT loves Healthy Eating – both signs, perhaps, of a preference for eating in.
• And what does each university use to quench their thirst? Harvard students like to drink wine (classy!), while Stanford students prefer beer (the social drink of choice).
• What about sports teams? MIT and Caltech couldn’t care less, though Harvard follows the Boston Red Sox, Stanford follows the San Francisco Giants (as well as their own Stanford Football and Stanford Cardinal), and Berkeley follows the Los Angeles Lakers (and the California Golden Bears).
• For sports themselves, MIT students like skiing; Stanford students like general exercise, fitness, and sports; Harvard students like baseball, ultimate frisbee, and skiing; and Berkeley students like snowboarding. Caltech, in a league of its own, enjoys table tennis and chess.
• What does each school think of social? Caltech students look for Social Advice. Berkeley students are interested in Social Media and Social Media Marketing. MIT, on the more technical side, wants Social Search. Stanford students, predictably, love the whole spectrum of social offerings, from Social Dance and The Social Network, to Social Psychology and Social Advice. (Interestingly, Caltech and Stanford are both interested in Social Advice, though I wonder if it’s for slightly different reasons.)
• What’s each school’s relationship with computers? Caltech students are interested in Computer Science, MIT hackers are interested in Computer Security, and Stanford students are interested in Human-Computer Interaction.
• Digging into the MIT vs. Caltech divide a little, we see that Caltech students really are more interested in the pure sciences (Physics, Science, Biology, Quantum Mechanics, Mathematics, Chemistry, etc.), while MIT students are more on the applied and engineering sides (Mechanical Engineering, Engineering, Distributed Databases, Cryptography, Computer Security, Biotechnology, Operating Systems, Compilers, etc.).
• Regarding programming languages, Caltech students love Haskell (hardcore purity!), while MIT students love Lisp.
• What does each school like to read, both offline and online? Caltech loves science fiction, xkcd, and The Onion; MIT likes Hacker News; Harvard loves journals, newspapers, and magazines (Huffington Post, the New York Times, Fortune, Wall Street Journal, the New Yorker, the Economist, and so on); and Stanford likes TechMeme.
• What movies and television shows does each school like to watch? Caltech likes Star Trek, the Colbert Report, and Futurama. MIT likes Fight Club (I don’t know what this has to do with MIT, though I will note that on my first day as a freshman in a new dorm, Fight Club was precisely the movie we all went to a lecture hall to see). Stanford likes The Social Network and Inception. Harvard, rather fittingly, likes Mad Men and Ted Talks.
• Let’s look at the startups each school follows. MIT, of course, likes Ksplice. Berkeley likes Yelp and Groupon. Stanford likes just about every startup under the sun (Instagram, Flipboard, Tumblr, Path, Color Labs, etc.). And Harvard, that bastion of hard-won influence and prestige? To the surprise of precisely no one, Harvard enjoys Klout.

Let’s end with a summarized view of each school:

• Caltech is very much into the sciences (Physics, Biology, Quantum Mechanics, Mathematics, etc.), as well as many pretty nerdy topics (Star Trek, Science Fiction, xkcd, Futurama, Starcraft II, etc.).
• MIT is dominated by everything engineering and tech.
• Stanford loves relationships (interpersonal relationships, people skills, love, network effects, sex, etiquette, dating and relationships, romance), health and appearance (fashion, fitness, nutrition, happiness), and startups (Instagram, Flipboard, Path, Color Labs, etc.).
• Berkeley, sadly, is perhaps too large and diverse for an overall characterization.
• Harvard students are fascinated by famous figures (Jimmy Fallon, Oprah Winfrey, Invaka Trump, Dalai Lama, David Lynch, Al Gore, Bill Gates, Barack Obama), and by prestigious newspapers, journals, and magazines (Fortune, the New York Times, the Wall Street Journal, the Economist, and so on). Other very fitting interests include Kiva, classical music, and Coldplay.

*I pulled about 400 followers from each school, and added a couple filters, to try to ensure that followers were actual attendees of the schools rather than general people simply interested in them. Topics are sorted using a naive Bayes score and filtered to have at least 5 counts. Also, a word of warning: my dataset was fairly small and users on Quora are almost certainly not representative of their schools as a whole (though I tried to be rigorous with what I had).

tl;dr See this movie visualization for a case study on how a post propagates through Quora.

How does information spread through a network? Much of Quora’s appeal, after all, lies in its social graph – and when you’ve got a network of users, all broadcasting their activities to their neighbors, information can cascade in multiple ways. How do these social designs affect which users see what?

Think, for example, of what happens when your kid learns a new slang word at school. He doesn’t confine his use of the word to McKinley Elementary’s particular boundaries, between the times of 9-3pm – he introduces it to his friends from other schools at soccer practice as well. A couple months later, he even says it at home for the first time; you like the word so much, you then start using it at work. Eventually, Justin Bieber uses the word in a song, at which point the word’s popularity really starts to explode.

So how does information propagate through a social network? What types of people does an answer on Quora reach, and how does it reach them? (Do users discover new answers individually, or are hubs of connectors more key?) How does the activity of a post on Quora rise and fall? (Submissions on other sites have limited lifetimes, fading into obscurity soon after an initial spike; how does that change when users are connected and every upvote can revive a post for someone else’s eyes?)

(I looked at Quora since I had some data from there already available, but I hope the lessons should be fairly applicable in general, to other social networks like Facebook, Twitter, and LinkedIn as well.)

To give an initial answer to some of these questions, I dug into one of my more popular posts, on a layman’s introduction to random forests.

Users, Topics

Before looking deeper into the voting dynamics of the post, let’s first get some background on what kinds of users the answer reached.

Here’s a graph of the topics that question upvoters follow. (Each node is a topic, and every time upvoter X follows both topics A and B, I add an edge between A and B.)

We can see from the graph that upvoters tend to be interested in three kinds of topics:

• Machine learning and other technical matters (the green cluster): Classification, Data Mining, Big Data, Information Retrieval, Analytics, Probability, Support Vector Machines, R, Data Science, …
• Startups/Silicon Valley (the red cluster): Facebook, Lean Startups, Investing, Seed Funding, Angel Investing, Technology Trends, Product Managment, Silicon Valley Mergers and Acquisitions, Asana, Social Games, Quora, Mark Zuckerberg, User Experience, Founders and Entrepreneurs, …
• General Intellectual Topics (the purple cluster): TED, Science, Book Recommendations, Philosophy, Politics, Self-Improvement, Travel, Life Hacks, …

Also, here’s the network of the upvoters themselves (there’s an edge between users A and B if A follows B):

We can see three main clusters of users:

• A large group in green centered around a lot of power users and Quora employees.
• A machine learning group of folks in orange centered around people like Oliver Grisel, Christian Langreiter, and Joseph Turian.
• A group of people following me, in purple.
• Plus some smaller clusters in blue and yellow. (There were also a bunch of isolated users, connected to no one, that I filtered out of the picture.)

Digging into how these topic and user graphs are related:

• The orange cluster of users is more heavily into machine learning: 79% of users in that cluster follow more green topics (machine learning and technical topics) than red and purple topics (startups and general intellectual matters).
• The green cluster of users is reversed: 77% of users follow more of the red and purple clusters of topics (on startups and general intellectual matters) than machine learning and technical topics.

More interestingly, though, we can ask: how do the connections between upvoters relate to the way the post spread?

Social Voting Dynamics

So let’s take a look. Here’s a visualization I made of upvotes on my answer across time (click here for a larger view).

To represent the social dynamics of these upvotes, I drew an edge from user A to user B if user A transmitted the post to user B through an upvote. (Specifically, I drew an edge from Alice to Bob if Bob follows Alice and Bob’s upvote appeared within five days of Alice’s upvote; this is meant to simulate the idea that Alice was the key intermediary between my post and Bob.)

Also,

• Green nodes are users with at least one upvote edge.
• Blue nodes are users who follow at least one of the topics the post is categorized under (i.e., users who probably discovered the answer by themselves).
• Red nodes are users with no connections and who do not follow any of the post’s topics (i.e, users whose path to the post remain mysterious).
• Users increase in size when they produce more connections.

Here’s a play-by-play of the video:

• On Feb 14 (the day I wrote the answer), there’s a flurry of activity.
• A couple of days later, Tracy Chou gives an upvote, leading to another spike in activity.
• Then all’s quiet until… bam! Alex Kamil leads to a surge of upvotes, and his upvote finds Ludi Rehak, who starts a small surge of her own. They’re quickly followed by Christian Langreiter, who starts a small revolution among a bunch of machine learning folks a couple days later.
• Then all is pretty calm again, until a couple months later when… bam! Aditya Sengupta brings in a smashing of his own followers, and his upvote makes its way to Marc Bodnick, who sets off a veritable storm of activity.

(Already we can see some relationships between the graph of user connections and the way the post propagated. Many of the users from the orange cluster, for example, come from Alex Kamil and Christian Langreiter’s upvotes, and many of the users from the green cluster come from Aditya Sengupta and Marc Bodnick’s upvotes. What’s interesting, though, is, why didn’t the cluster of green users appear all at once, like the orange cluster did? People like Kah Seng Tay, Tracy Chou, Venkatesh Rao, and Chad Little upvoted the answer pretty early on, but it wasn’t until Aditya Sengupta’s upvote a couple months later that people like Marc Bodnick, Edmond Lau, and many of the other green users (who do indeed follow that first set of folks) discovered the answer. Did the post simply get lost in users’ feeds the first time around? Was the post perhaps ignored until it received enough upvotes to be considered worth reading? Are some users’ upvotes just trusted more than others’?)

For another view of the upvote dynamics, here’s a static visualization, where we can again easily see the clusters of activity:

Fin

There are still many questions it would be interesting to look at; for example,

• What differentiates users who sparked spikes of activity from users who didn’t? I don’t believe it’s simply number of followers, as many well-connected upvoters did not lead to cascades of shares. Does authority matter?
• How far can a post reach? Clearly, the post reached people more than one degree of separation away from me (where one degree of separation is a follower); what does the distribution of degrees look like? Is there any relationship between degree of separation and time of upvote?
• What can we say about the people who started following me after reading my answer? Are they fewer degrees of separation away? Are they more interested in machine learning? Have they upvoted any of my answers before? (Perhaps there’s a certain “threshold” of interestingness people need to overflow before they’re considered acceptable followees.)

But to summarize a bit what we’ve seen so far, here are some statistics on the role the social graph played in spreading the post:

• There are 5 clusters of activity after the initial post, sparked both by power users and less-connected folks. In an interesting cascade of information, some of these sparks led to further spikes in activity as well (as when Aditya Sengupta’s upvote found its way to Marc Bodnick, who set off even more activity).
• 35% of users made their way to my answer because of someone else’s upvote.
• Through these connections, the post reached a fair variety of users: 32% of upvoters don’t even follow any of the post’s topics.
• 77% of upvotes came from users over two weeks after my answer appeared.
• If we look only at the upvoters who follow at least one of the post’s topics, 33% didn’t see my answer until someone else showed it to them. In other words, a full one-third of people who presumably would have been interested in my post anyways only found it because of their social network.

So it looks like the social graph played quite a large part in the post’s propagation, and I’ll end with a big shoutout to Stormy Shippy, who provided an awesome set of scripts I used to collect a lot of this data.

Introduction

Suppose you have the following set of sentences:

• I like to eat broccoli and bananas.
• I ate a banana and spinach smoothie for breakfast.
• Chinchillas and kittens are cute.
• My sister adopted a kitten yesterday.
• Look at this cute hamster munching on a piece of broccoli.

What is latent Dirichlet allocation? It’s a way of automatically discovering topics that these sentences contain. For example, given these sentences and asked for 2 topics, LDA might produce something like

• Sentences 1 and 2: 100% Topic A
• Sentences 3 and 4: 100% Topic B
• Sentence 5: 60% Topic A, 40% Topic B
• Topic A: 30% broccoli, 15% bananas, 10% breakfast, 10% munching, … (at which point, you could interpret topic A to be about food)
• Topic B: 20% chinchillas, 20% kittens, 20% cute, 15% hamster, … (at which point, you could interpret topic B to be about cute animals)

The question, of course, is: how does LDA perform this discovery?

LDA Model

In more detail, LDA represents documents as mixtures of topics that spit out words with certain probabilities. It assumes that documents are produced in the following fashion: when writing each document, you

• Decide on the number of words N the document will have (say, according to a Poisson distribution).
• Choose a topic mixture for the document (according to a Dirichlet distribution over a fixed set of K topics). For example, assuming that we have the two food and cute animal topics above, you might choose the document to consist of 1/3 food and 2/3 cute animals.
• Generate each word w_i in the document by:
  • First picking a topic (according to the multinomial distribution that you sampled above; for example, you might pick the food topic with 1/3 probability and the cute animals topic with 2/3 probability).
  • Using the topic to generate the word itself (according to the topic’s multinomial distribution). For example, if we selected the food topic, we might generate the word “broccoli” with 30% probability, “bananas” with 15% probability, and so on.

Assuming this generative model for a collection of documents, LDA then tries to backtrack from the documents to find a set of topics that are likely to have generated the collection.

Example

Let’s make an example. According to the above process, when generating some particular document D, you might

• Pick 5 to be the number of words in D.
• Decide that D will be 1/2 about food and 1/2 about cute animals.
• Pick the first word to come from the food topic, which then gives you the word “broccoli”.
• Pick the second word to come from the cute animals topic, which gives you “panda”.
• Pick the third word to come from the cute animals topic, giving you “adorable”.
• Pick the fourth word to come from the food topic, giving you “cherries”.
• Pick the fifth word to come from the food topic, giving you “eating”.

So the document generated under the LDA model will be “broccoli panda adorable cherries eating” (note that LDA is a bag-of-words model).

Learning

So now suppose you have a set of documents. You’ve chosen some fixed number of K topics to discover, and want to use LDA to learn the topic representation of each document and the words associated to each topic. How do you do this? One way (known as collapsed Gibbs sampling) is the following:

• Go through each document, and randomly assign each word in the document to one of the K topics.
• Notice that this random assignment already gives you both topic representations of all the documents and word distributions of all the topics (albeit not very good ones).
• So to improve on them, for each document d…
  • Go through each word w in d…
    • And for each topic t, compute two things: 1) p(topic t | document d) = the proportion of words in document d that are currently assigned to topic t, and 2) p(word w | topic t) = the proportion of assignments to topic t over all documents that come from this word w. Reassign w a new topic, where we choose topic t with probability p(topic t | document d) * p(word w | topic t) (according to our generative model, this is essentially the probability that topic t generated word w, so it makes sense that we resample the current word’s topic with this probability). (Also, I’m glossing over a couple of things here, in particular the use of priors/pseudocounts in these probabilities.)
    • In other words, in this step, we’re assuming that all topic assignments except for the current word in question are correct, and then updating the assignment of the current word using our model of how documents are generated.
• After repeating the previous step a large number of times, you’ll eventually reach a roughly steady state where your assignments are pretty good. So use these assignments to estimate the topic mixtures of each document (by counting the proportion of words assigned to each topic within that document) and the words associated to each topic (by counting the proportion of words assigned to each topic overall).

Layman’s Explanation

In case the discussion above was a little eye-glazing, here’s another way to look at LDA in a different domain.

Suppose you’ve just moved to a new city. You’re a hipster and an anime fan, so you want to know where the other hipsters and anime geeks tend to hang out. Of course, as a hipster, you know you can’t just ask, so what do you do?

Here’s the scenario: you scope out a bunch of different establishments (documents) across town, making note of the people (words) hanging out in each of them (e.g., Alice hangs out at the mall and at the park, Bob hangs out at the movie theater and the park, and so on). Crucially, you don’t know the typical interest groups (topics) of each establishment, nor do you know the different interests of each person.

So you pick some number K of categories to learn (i.e., you want to learn the K most important kinds of categories people fall into), and start by making a guess as to why you see people where you do. For example, you initially guess that Alice is at the mall because people with interests in X like to hang out there; when you see her at the park, you guess it’s because her friends with interests in Y like to hang out there; when you see Bob at the movie theater, you randomly guess it’s because the Z people in this city really like to watch movies; and so on.

Of course, your random guesses are very likely to be incorrect (they’re random guesses, after all!), so you want to improve on them. One way of doing so is to:

• Pick a place and a person (e.g., Alice at the mall).
• Why is Alice likely to be at the mall? Probably because other people at the mall with the same interests sent her a message telling her to come.
• In other words, the more people with interests in X there are at the mall and the stronger Alice is associated with interest X (at all the other places she goes to), the more likely it is that Alice is at the mall because of interest X.
• So make a new guess as to why Alice is at the mall, choosing an interest with some probability according to how likely you think it is.

Go through each place and person over and over again. Your guesses keep getting better and better (after all, if you notice that lots of geeks hang out at the bookstore, and you suspect that Alice is pretty geeky herself, then it’s a good bet that Alice is at the bookstore because her geek friends told her to go there; and now that you have a better idea of why Alice is probably at the bookstore, you can use this knowledge in turn to improve your guesses as to why everyone else is where they are), and eventually you can stop updating. Then take a snapshot (or multiple snapshots) of your guesses, and use it to get all the information you want:

• For each category, you can count the people assigned to that category to figure out what people have this particular interest. By looking at the people themselves, you can interpret the category as well (e.g., if category X contains lots of tall people wearing jerseys and carrying around basketballs, you might interpret X as the “basketball players” group).
• For each place P and interest category C, you can compute the proportions of people at P because of C (under the current set of assignments), and these give you a representation of P. For example, you might learn that the people who hang out at Barnes & Noble consist of 10% hipsters, 50% anime fans, 10% jocks, and 30% college students.

Real-World Example

Finally, I applied LDA to a set of Sarah Palin’s emails a little while ago (see here for the blog post, or here for an app that allows you to browse through the emails by the LDA-learned categories), so let’s give a brief recap. Here are some of the topics that the algorithm learned:

• Trig/Family/Inspiration: family, web, mail, god, son, from, congratulations, children, life, child, down, trig, baby, birth, love, you, syndrome, very, special, bless, old, husband, years, thank, best, …
• Wildlife/BP Corrosion: game, fish, moose, wildlife, hunting, bears, polar, bear, subsistence, management, area, board, hunt, wolves, control, department, year, use, wolf, habitat, hunters, caribou, program, denby, fishing, …
• Energy/Fuel/Oil/Mining: energy, fuel, costs, oil, alaskans, prices, cost, nome, now, high, being, home, public, power, mine, crisis, price, resource, need, community, fairbanks, rebate, use, mining, villages, …
• Gas: gas, oil, pipeline, agia, project, natural, north, producers, companies, tax, company, energy, development, slope, production, resources, line, gasline, transcanada, said, billion, plan, administration, million, industry, …
• Education/Waste: school, waste, education, students, schools, million, read, email, market, policy, student, year, high, news, states, program, first, report, business, management, bulletin, information, reports, 2008, quarter, …
• Presidential Campaign/Elections: mail, web, from, thank, you, box, mccain, sarah, very, good, great, john, hope, president, sincerely, wasilla, work, keep, make, add, family, republican, support, doing, p.o, …

Here’s an example of an email which fell 99% into the Trig/Family/Inspiration category (particularly representative words are highlighted in blue):

And here’s an excerpt from an email which fell 10% into the Presidential Campaign/Election category (in red) and 90% into the Wildlife/BP Corrosion category (in green):