Python Guy

Probabilities December 5, 2011

I used to play a game called Hamurabi, a very ancient game.

The way it works is you have a population, an amount of land, and an amount of grain. Every year, you need to choose how much grain to plant, how much grain to feed the people, and how many acres of land to buy or sell.

Each year, a certain amount of grain would grow, people would starve if they didn’t get enough food, and rats would eat your surplus grain. Random events would cause disasters to make things interesting.

The game was remarkably simple but surprisingly addictive.

I figured it was time to write a new version for the modern age, when memory and CPU are no longer real constraints to most of the problems we faced programming 30 years ago. I chose the Hamurabi game as inspiration.

Of course, a real simulation of a primitive society could track each individual within the population. People would age, and then die from starvation, random events, or just old age.

Using the US Census data, I found some actuarial tables that gave a percentage chance someone would die at a particular age.

I stored the number of people in each age group, separated by sex, in a simple array. The index was the age.

Now all I needed to do was find a scalable way of calculating how many people die in each age group every year.

The simplest way to do this is to roll the dice for each person each year. Simply use the random modules’ random() function to get a number between 0.0 and 1.0. If this is left than the probability they will die, then they die.

Of course, if you have a thousand people, you have to call random() a thousand times. A million people need it called a million times. This would slow the program down as your population grew.

I vaguely recalled some math from my college years that should help me work these things out more quickly. Reading some of my old textbooks, I discovered that that Binomial Distribution gives you the exact probability of seeing x events in population of n where each has a probability of p.

However, the binomial distribution relies on combinatorials, which rely on factorials, which are notoriously difficult to calculate precisely for large numbers.

Luckily, this is not a new problem. Approximating the Binomial Distribution leads to the Normal or Gauss Distribution. This relies on simple exponentials and Python’s random module already has a gauss() and a normalvariate() built in.

What do you use for mu and sigma, the mean and standard deviation? The mean is simply n*p, and the standard deviation is the square root of n*p*(1-p).

The normal distribution really doesn’t give you good approximations when n or n*p is small. If n is really small (there are only a handful of people), I can just roll random() that many times. But if n is large, and p is very small (generally, for younger ages), then my approximation breaks down.

For this, the Poisson Distribution is useful. The Poisson Distribution can give you a very good approximation of the Binomial Distribution for small numbers of deaths. What I do is roll the dice with random(), and walk from 0 deaths upward, subtracting the probability that that many deaths would occur given the Poisson Distribution. When I’ve exceeded the roll, then that’s the number of deaths that occur.

Using these three methods of calculating deaths lead to a very scalable algorithm for calculating the number of deaths in a population of 1, 10, hundreds of even millions of people. I ran my simulation for several hundred years until I saw hundreds of billions of people, with tens of billions dying each year.

I think this shows why I like Python. Never once did I have to think very hard about how to express my ideas in Python. I was completely free to consider the algorithm. Most of my time was spent reading math books and translating complicated formulas into simple algorithms, and testing that it did what I expected.