Wednesday, October 22nd, 2008

A nifty little election time dynamic program

Take a look at xkcd's election predictor. It scrapes a bunch of election outcome probabilities per state from intrade and uses them to provide a prediction for the overall election. It does so by using a Monte Carlo simulation, given the probabilities, it runs a mock election assuming each state outcome is independent a whole bunch of times with those probabilities and see what happens. In his code, he runs the simulation 100,000 times. On my laptop, this takes about 4 seconds and yields only an approximation with error of about .001. Aww, poor physicist, if only Randall was a computer scientist, maybe he'd see how to compute it exactly and much more quickly.

So, how can we compute this efficiently? The key to my program running quickly is that total number of electoral votes (538) is much smaller than the total number of possible election outcomes (2^50).

Imagine we have a partially computed an election in which k states have already voted, and that pk,i is the probability that the democrats have won exactly i electoral votes after the first k elections have run. Furthermore, let vk be the number of electoral votes for state k and ek be the probability state that state k elects a democrat.

How can we compute pk,i? Well, this simple little recurrence will do it.
pk,i = ek * pk-1,i - vk + (1 - ek) * pk-1, i. Basically, with probability ek the state votes democrat and has i - vk votes left, otherwise, with probability 1 - ek it votes republican (assuming two party system), and the democrat still has i votes left. For this problem, we only need about 539 * 51 table entries (we could even do something clever and push the needed space down to 539 entries, but I'll leave that as an exercise for the reader).

After I implemented this, the computation runs much more quickly than one second and provides an exact answer.

Looking at the code a bit more, I found a (slightly refactored) snippet of code like this. I assume the code is supposed to output the outcomes in order. Take a look at the comment, is it simpler than sorting? Is it even correct? I challenge both accounts.

def render_outcome(pdem, prep, ptie):
    dstring="Obama:  <span style=\"color: #0000FF\">"+str(round(pdem,1))+"</span>"
    rstring="McCain: <span style=\"color: #FF0000\">"+str(round(prep, 1))+"</span>"
    tstring="Tie: <span style=\"color: #888888\">"+str(round(ptie, 1))+"</span>"

    if pdem>prep and prep>ptie:
        return dstring+"\n"+rstring+"\n"+tstring
    if prep>pdem and pdem>ptie:
        return rstring+"\n"+dstring+"\n"+tstring
    #just in case ... (easier to do this than a sort)
    if pdem>ptie and ptie>prep:
        return dstring+"\n"+tstring+"\n"+rstring
    if prep>ptie and ptie>pdem:
        return rstring+"\n"+tstring+"\n"+dstring
    if ptie>pdem and pdem>prep:
        return tstring+"\n"+dstring+"\n"+rstring
    if ptie>prep and prep>pdem:
        return tstring+"\n"+rstring+"\n"+dstring
    return tstring+"\n"+rstring+"\n"+dstring


Here is my alternative implementation.

def render_outcome2(pdem, prep, ptie):
    dstring="obama:  <span style=\"color: #0000ff\">"+str(round(pdem,1))+"</span>"
    rstring="mccain: <span style=\"color: #ff0000\">"+str(round(prep, 1))+"</span>"
    tstring="tie: <span style=\"color: #888888\">"+str(round(ptie, 1))+"</span>"
    l = [(pdem, dstring), (prep, rstring), (ptie, tstring)]
    return '\n'.join(pair[1] for pair in sorted(l, reverse=True))


Mine definitely seems simpler. It relies on the natural sorting order of python tuples to get the messages sorted in the right order.

Is his implementation correct? Well.. notice all of those < operators (not <=). What happens with ties?

>>> print states.render_outcome(.4, .3, .3)
Tie: 0.3
McCain: 0.3
Obama:  0.4


Uh oh..

In all fairness, quoting the page, Randall says "I made this tool to help me understand the race, especially on election night." I am sure he just wanted to get things done, and not have some nerd nitpick at all of the code. The Monte Carlo simulation is a bit easier to code than the dynamic program I posted and it gets things done. His code basically works. I don't think he actually sucks at programming, I just wanted to put some blood on the pages for reddit. Furthermore, I was thinking about using this problem as an interview question, but after trying it on a few of my coworkers (who all have at least a BS in computer science from a nice university), I think it's a bit too hard.
(12 comments | Leave a comment)

Friday, March 21st, 2008

Travelling salesman, xkcd, and a month old reddit post.

From a reddit post I made a month ago.


There are n! paths in a graph of size n. I can find the optimal TSP solution in O(n^2 2^n) time. Therefore, it is possible to do better than examine every path to find a TSP solution.

http://www.algorithmist.com/index.php/Traveling_Salesperson_Problem


And today's xkcd comic.



Suspicious?

I just wish I had a modicum of artistic ability, I swear I've got about five ideas for good xkcd style comics.
(2 comments | Leave a comment)

Sunday, January 6th, 2008

The right way to understand repeated squaring

So, I wish I could say I was over this, but I am still having bouts of sadness and self-induced frustration. But at least I am having dreams about algorithms, which you know, kind of rocks. Even if my ability to reason in my dreams is laughably bad.

In my most recent dream, I was arguing in a classroom about the "right" way to understand repeated squaring. This was the argument that I wanted to make, but failed to express clearly in my dream.

I wrote most of the algorithmist article on repeated squaring awhile ago. The way I used to understand repeated squaring is explained by example in the write up. It basically comes down to writing the exponent in binary, accumulating the appropriate terms for each position the binary expansion of the exponent by repeatedly squaring the last position, and then finally multiplying out all the terms where there is a bit set. This works and provides some insight into how the algorithm works.

However, this is the wrong way to understand the algorithm.

The right way is expressed beautifully by this equation.

x^y = (x^{y~{\rm div}~2})^2 \cdot x^{y\bmod 2}

What does the first part equation say?

(x^{y~{\rm div}~2})

Well, it says divide the problem into problem of approximately half the size. This smaller problem can be solved recursively. The division here is the C-style truncating integer division. The square then uses the information from the subproblem to compute the larger problem. The last part of the equation,

x^{y\bmod 2}

Simply handles the case that the division truncated exponent and we would have otherwise lost a multiple of the base.

On a technical note, from either view of the algorithm, it's pretty easy to
see that the number of total number of multiplications required is no worse than 2lg(y).

What is the key difference between the two views of the algorithm? The first is iterative, understanding the computation involves understanding every minute detail of the accumulation loop and it's relation to the binary expansion of the exponent. The second view uses recursion to leverage abstraction. Instead of understanding the whole process, you simply need to understand how to make the problem smaller and how to put the pieces together on a high level. There is really a parallel here between iterative and functional programming. Ironically, I am coming out on the functional programming side.

So why is the second view better than the first view? The higher level, more abstracted view is much easier to generalize. Instead of computing x^y efficiently, let's consider the problem of computing a^b * x^y. Let n = max(b, y). This can be done straight-forwardly with two applications of the repeated squaring algorithm in 4lg(n) + some constant k multiplications. However, a modified algorithm can solve the problem with 2lg(n) + k multiplications. Credit for this problem goes to my sister's "Reasoning about Computation" class. It's the best class I never took. Sigh, sometimes I wish I did CS at Princeton.

How do you do it?

Here is some empty space in case you want to try to figure it out yourself. Hint: Try to use divide and conquer.















Well, to compute a^b * x^y, first compute a^(b/2) * x^(y/2) recursively. Now square it. There are just 4 cases depending on the truncation of b and y. If b and y were both truncated, we need an additional multiple of a * x, which can be precomputed ahead of time for cost of a single multiplication at this recursive step. Otherwise, if only b was truncated, we just need an additional multiple of a. Likewise, if only y was truncated, we need an additional multiple of x. Otherwise, no additional multiplications are needed. Simple, elegant, clever. Gotta love it.

Exercise for the readers: Further generalize the algorithm. Show how to compute a_1 ^ b_1 * a_2 ^ b_2 * ... a_j ^ b_j in 2 lg(n) + 2^j + k multiplications, where n = max(b_1, b_2, ... b_j) and k is a constant.

I spent quite a bit of time trying to solve this problem, mostly stuck because I was trying to use the first understanding of the problem and I couldn't see how to leverage the binary expansion of both exponents simultaneously.
(5 comments | Leave a comment)