Probability Problems Programmatically

Birthday Problem

How big of a group of people do you need so that there is a 50% chance 2 of them share a birthday? How about 99.9%? Of course, 100% is 366 people by the pigeon principle, but the small numbers for lower probabilities is surprising.

The key is that birthday comparisons are made between every pair of individuals and not just between one individual and everyone else. Then the solution becomes p(n) = 1 - ₃₆₅P_n/365ⁿ. At n=23, the probability of two individuals having the same birthday is about 50.7%.

If we didn't know how to arrive at this solution, we could use a simple Monte Carlo simulation to estimate the probability:

Monte Hall Problem
You are on a game show and are given the opportunity to select one of 3 doors. Behind one door, known to the host, is a major prize. Behind the other 2 doors are goats. Once you've made a selection the host will open a door showing a goat. Now you have the opportunity to select a different door. Should you switch or keep your original guess?

This problem is very counter-intuitive and seems to cause disagreement among otherwise intelligent people even after the correct solution is known. By switching doors, you dramatically increase your chance of winning the prize from 1/3 to 2/3.

Again, instead of worrying too much about how to arrive at the correct solution we can simulate it:

Determining if a coin is fair

Out of 1000 coin flips, how many heads can you get and still assume with 95% confidence that it's a fair coin?

Coins flips are a binomial distribution, but can be approximated as a normal distribution based on the central limit theorem. A 95% confidence interval is ± 1.96 standard deviations. The variance in 1000 flips is 15.8. 530 heads is (530-500)/15.8 or ± 1.9 standard deviations.

Rain in Seattle Problem

You call 3 friends and ask them if it's raining. Each friend has a 1/3 chance of lying and 2/3 chance of telling the truth. All 3 friends say it's raining. What is the probability it's actually raining?
If we say it's raining 10% of the time, then our probability is .1 * (2/3)3 / (.1 * (2/3)3 + .9 * (1/3)3) by Baye's theorem. The answer is 47%.

These simple Monte Carlo simulations are a way to get an approximate answer when we don't know how to (or don't want to) calculate the exact answer. These were all just a few lines of code, took a few minutes to write, and run nearly instantly. It's a great way to verify the above answers as well. And to get a more exact answer, you just trade off the quick run time with doing more than 10,000 simulations.

Linear Programming

Linear programming (LP) is a technique for the optimization of a linear objective function. Integer programming (IP) has the additional constraint of the unknown variables being integers. I vaguely remember these from school, but they are terms I see every so often so I wanted to code something up and get reacquainted.

Many LP solvers exist and a quick search for a Python library turned me onto PuLP. They already had a Sudoku example, which is what I wanted to do, so I've turned that into a Jupyter notebook (getting one of these embedded in a post is something else I've been wanting to try) to make it a bit easier to follow.

LP has polynomial time solutions (see Karmarkar's algorithm), while IP is NP-hard (see branch-and-bound). With Sudoku being a relatively small problem it doesn't seem to matter. These puzzles are easily solvable in under 1 second.

The linear programming I remember was graphing several lines and getting a solution where they intersect. Applying the technique to Sudoku is way more fun. And there are obviously other more useful applications like portfolio optimization, airline crew scheduling, or vehicle routing.

While working on this I also came across Google Optimization Tools which I haven't heard of before. I don't know how popular PuLP is, so it could be another, perhaps better supported, alternative.

UPDATE:

CVXPY also looks promising as demonstrated in Optimization with Python: How to make the most amount of money with the least amount of risk.

Integer vs Linear Programming in Python has a nice comparison of ML to linear programming that will help you decide when to use each.

Python 3 asyncio

A while back I wrote a few posts about asynchronous programming:

• Asynchronous Non-blocking I/O Java Echo Server
• Fixing My Asynchronous Non-blocking I/O Callback Hell With Monads
• Scala Asynchronous IO Echo Server
• Reactive Actors

So when I learned that Python 3.5 has added async and await operations I knew I had to check it out to see how it compares. asyncio describes itself as "infrastructure for writing single-threaded concurrent code using coroutines, multiplexing I/O access over sockets and other resources, running network clients and servers, and other related primitives".

Coroutines in Python are similar to generators, but coroutines (a function definition using async def) can control where execution continues after the yield (replaced by await). You await a another coroutine or a future. The coroutine approach can replace callbacks.

If we check the time it takes for the program to run, it's 3 seconds (the longest slow operation) and not 6 seconds (the sum of all the slow operations).

Slow operation sleep 1 complete

Slow operation sleep 2 complete

Slow operation sleep 3 complete

Completed in 3.02 seconds

It's similar if we want to get results from each task. We're even able to get them as they are available. 

Slow operation sleep 1 complete

Got result 1

Slow operation sleep 2 complete

Got result 2

Slow operation sleep 3 complete

Got result 3

Completed in 3.00 seconds

With Python's GIL you've never really been able to run multiple threads in parallel. You've had to run concurrent code in multiple processes to leverage multiple CPU cores. With asyncio we are able to at least make single process IO-bound tasks execute faster because it switches between tasks to bypass GIL contention.

For CPU-bound code you still have to use multiple processes to parallelize your code. Python's parallel API limits your ability to use results from these tasks as they become available, as we did in the example above. Waiting for a process to finish and getting a future result both block. asyncio gives us a way to unify the concurrent and parallel APIs.

Slow operation sum 10000000 complete

Got result 49999995000000

Slow operation sum 20000000 complete

Got result 199999990000000

Slow operation sum 30000000 complete

Got result 449999985000000

Completed in 2.91 seconds

There is some overhead in creating the processes but a single 3 second task is only slightly faster than parallel 1, 2, and 3 second tasks.

Slow operation sum 30000000 complete

Got result 449999985000000

Completed in 2.66 seconds

This hints at another use. Given that asyncio uses a single thread, if there are too many IO-bound tasks or if any of them consumes too much CPU, we can overwhelm it. For such a situation it is possible to create multiple processes, each with its own asyncio event loop.

1 process and 3 tasks:

Got 3 results in 6.20 seconds

1 process and 15 tasks:

Got 15 results in 23.94 seconds

The simulated mix of IO-bound and CPU-bound code is interesting. It's slower than just the CPU-bound code on it's own, but faster than the sum of the IO-bound and CPU-bound code. We see some asyncio benefits up to around 15 tasks and then it levels off.

1 process and 30 tasks:

Got 30 results in 46.58 seconds

4 processes and 15*4 tasks:

Got 60 results in 34.11 seconds

8 processes and 15*8 tasks:

Got 120 results in 52.81 seconds

We similarly see some benefit of parallelizing the code up the number of CPU cores I have, then the time is increasing linearly as we would expect.

16 processes and 15*16 tasks:

Got 120 results in 105.98 seconds

If you are looking for more details the comprehensive Python Concurrency series of articles has additional examples like this and goes into more detail on some of the underlying concepts.