8 Temmuz 2012 Pazar

Counting and Number Systems

To contact us Click HERE
We all learned how to count as young children: "one, two, three". Iflike me you watched a lot of Sesame Street when you were small, thenyou also learned the Spanish words for those same concepts: "uno, dos,tres". And that pretty much sums up counting, right?

Well, not quite. Today, I would like to try something a littledifferent. Instead of discussing a specific application of math tobusiness, science, or art, I want to go back to the very beginning andsay a few words about "counting" that illustrate the relationshipbetween mathematical abstractions and the real world.

How did early humans first develop mathematics? How do/should childrenlearn about math? Certainly one concept at the heart of mathematicalthinking is counting. Hunters count animals, gatherers countfruits, children count marbles, merchants count inventory (both goodsand coins): pretty much everyone needs to count.

Pretend you are just learning how to count, and you have to make it upas you go along, without cultural knowledge to speed up theprocess. What do you do?

Perhaps you start by making marks in one-to-one correspondence withthe objects you are counting. Thus if you have one apple, you make onemark, as in Figure A. If you have two apples, you make two marks, asin Figure B. If you have six apples, you make six marks, as inFigure F.

This is a good start: it is verifiable (you can visually line up themarks and the apples and confirm that they are in a one-to-onecorrespondence), and it suggests the abstract concept of "number",because "five marks" could mean "five apples" or "five rabbits" orwhatever.

However, it is not that great as a system for organizing yourbusiness. It is a bit like making a map of your town at a scale of oneinch equals one inch: the map has to be as large as your town.What you really want is data compression: a map that fits inyour pocket so you have it when you go out walking, yet summarizesenough about the town to be useful in helping you find your way.

Can you visually distinguish between 5 and 6 marks (figures E and Fabove)? What about between 36 and 37 marks (not shown): short ofre-counting them, it is tough for the human visual system to processthat many marks: the Monotonous Repetition is Unsatisfyingeffect comes in, and our minds get bored.

So we need a better system. The first step we learned as childrenis to use a different symbol to represent groups of five. Forinstance, we can use every fifth "stick" to "tie up a bundle". Underthis system we can more easily distinguish numbers like 4, 5, 6, and 7:

However, it is still difficult to distinguish 101 from 106, because westill have monotonous repetition, just now it is bundles of five thatrepeat.

Nonetheless, we can generalize this approach. Every time we have 5 bundles(of five), let's replace them by yet another new symbol, which standsfor 25 sticks. Since drawing little stick figures is tedious, I amalso going to start writing them using letters rather than pictures,but that's just for my typing convenience. In fact, to make clear thatthese are still "pictures", I will use rather uncommon letters: S willmean a single "stick", B will mean a "bundle of five", and "G" will meana group of five bundles, i.e. 25 sticks. With this system, I can nowcount fairly high without making the symbols too monotonouslyrepetitive. Here are the symbols for a variety of differentnumbers. I have labeled them with their conventional English word names soyou can check your ability to decode them, but remember that at thispoint, we have no English words, no decimal system, only these GBSsymbols.

Symbolic PictographEnglish Interpretation
Sone
SStwo
SSSthree
SSSSfour
Bfive
BSsix
BSSseven
BBten
BBSeleven
BBBfifteen
BBBBtwenty
BBBBSSSStwenty four
Gtwenty five
GStwenty six
GBSthirty one
GGBBSSsixty two
GGGGSone hundred one

At this point, we have more or less replicated "Roman Numerals". Theactual Roman Numeral system is a bit more arcane than this, but sinceit is a mathematical dead-end, I am not going to bore you by walkingthrough it. The fatal flaw in both Roman Numerals and our little "GBS"system is that arithmetic is extremely difficult. We have succeeded ingiving (somewhat) short names to (small positive whole) numbers, in asystematic fashion. This makes it relatively easy to count. Butonce we can count, we want to do other things. Most urgently, we wantto add, so that if you have 27 apples and I have 32, we learnthat together we have accumulated a total of 59.

How would addition work in the GBS system? We could try to concatenatethe numbers for 27 and 32, namely "GSS" and "GBSS", which gives (aftersorting), "GGBSSSS", or 59. This is actually correct! Unfortunately,this simple rule fails in lots of other cases. For instance, if we add27 and 33, we get "GGBSSSSS", which is not even a valid number in oursystem. We see the problem: we need to "promote" the five S charactersinto a B, forming "GGBB" or 60.

What about adding 24 and 24?We get "BBBBBBBBSSSSSSSS", which is a huge monotonous mess. Now wehave to group five S characters into a B, forming "BBBBBBBBBSSS". Nowwe group five B characters into a G, forming "GBBBBSSS", or 48.Is this really how you want to do addition? Imagine tryingsubtraction, multiplication, or division!!

Also, if this is your system for writing numbers, how will you devise amechanical "adding machine" to help you? Notice that adding twonumbers sometimes creates a result that is shorter than either of thetwo inputs (e.g. add ten and fifteen). That is going to seriouslycomplicate any kind of system of gears and levers you might try to build.

So now we are in a position to better state what we want a numbersystem to achieve. Not only must it provide a systematic wayto count (i.e. it must give short names to small numbers, andbe able to systematically name any number), but it must alsosupport a simple, systematic way to do addition.It would be great if it also made it equally easy to multiply and divide, butthat may be asking too much - that's why people invented logarithmsand slide rules to supplement adding machines.

It turns out that there are many ways to accomplish thisgoal. Our familiar "decimal system" is one way, but since it is sofamiliar, I am going to present two others instead that are also quitehandy to know about in this age of computers: the "binary system" andthe "hexadecimal system". There are even more exotic systems. Forinstance, Knuth(The Art ofComputer Programming, Volume 2) discusses a perfectly reasonableQuater-imaginarybase number system.

The decimal, binary, hexadecimal and Quater-imaginary systemsall share a key feature: positional notation. This is anotherconceptual leap forward. Although I sorted the letters in my "GBS"system to put the "important" (big) letters first, that was purely atypographical convention which had no real impact: "GBS" and "SBG"both represent the number that we ordinarily call 31, because both arereally just short hand for drawing 31 sticks.

By contrast, in positional number systems, the order inwhich we write the symbols matters. This means that we can achievehigher amounts of data compression, because position can carryinformation implicitly, without us having to spell itout. However, data compression reduces redundancy.Unfortunately, this makes math hard for many people: less redundancymeans greater attention to detail is required to get it right. Aperson whose brain has trouble distinguishing between "GBS" and "GSB"is likely to have trouble with arithmetic.

Ordinary English text contains a great deal ofredundancy. Psychologists have done experiments where they show peoplethe top half of the letters in a sentence, or the bottom half, orevery third letter, and people can often figure out the meaning. If Iwrite "rss r rd, vlts r bl", you might be able to guess that I meant"roses are red, violets are blue". But in math, missing characters arefatal: each one carries its own crucial meaning - there is no analogyto "deleting vowels". There is no question that the numbers 5027 and507 are completely different, and both are contextuallyplausible, or equally "grammatical". You cannot drop (or swap) symbolsin math.

The binary system uses just two fundamental symbols,conventionally written 0 and 1. "Bi" means "two", like in"bicycle". In contrast, "deci" means ten and "hexideci"means 16, so those systems use more fundamental symbols ("digits").The decimal system uses the ten symbols "0123456789". Thehexadecimal system uses these and six more, namely the letters "abcdef".

Since binary uses just two symbols, it is simple to learn. Thesimplicity also facilitates building "adding machines" (nowadays knownas "computers").

Binary is a positional system. When we count, we start with 1. For thenext number, two, we do not have a '2' symbol, so we write 10instead. At this point you are probably saying, wait a minute, 10means "ten", not "two", and you feel hopelessly confused.

So, for purposes of explanation, let's not use 0 and 1 as the binarysymbols: let's use A and B. When we write a binary number such asBABA, the positions matter. The last (rightmost) symbolrepresents individual sticks/apples/rabbits, with A meaning zero and Bmeaning one. The second to last symbol represents bundles oftwo sticks: A still means zero, but B now means one bundle,i.e. two sticks. That's why "BA" stands for"two". The third to last symbol represents groupsof four sticks: A still means zero, but B now means one group,i.e. four sticks. The fourth to last symbol representsgroups of eight sticks, and so forth. Thus "BABA" means "ten", becauseit represents eight plus two.

Here is the same table of numbers we used to describe the GBS system,now converted to binary.

Binary SymbolEnglish Interpretation
Bone
BAtwo
BBthree
BAAfour
BABfive
BBAsix
BBBseven
BABAten
BABBeleven
BBBBfifteen
BABAAtwenty
BBAAAtwenty four
BBAABtwenty five
BBABAtwenty six
BBBBBthirty one
BBBBBAsixty two
BBAABABone hundred one

For example, the "one hundred one" case works like this.

There are many good things about this table: larger numbers (inmagnitude) are also longer (in binary form) than smaller ones, so thepattern is more systematic, and easier to mechanize, than for the GBSsystem. Moreover, it is easy to do addition: just apply the samealgorithm you learned in elementary school! Starting from theright-most column, you add, and if the result is too big to fit, you"carry" a "B" ("one") to the column to the left.

The nice part is you don't have to memorize such a big table ofaddition facts! School children have to memorize all possiblesingle-digit addition facts, such as "9+6 = 15". But in binary, theONLY non-trivial single-bit addition fact is "B+B = BA", i.e. the rulethat if both "bits" are "B", you write "A" for the sum and you "carry"the "B".

If you already had two B's in a column, and you received a third oneas a carry, you use the fact that
B+B+B = (B+B)+B = BA+B = BB,
so you would put B as the answerand carry another B to the left.

For instance, to add 10 and 6 to get 16, we work in binary as shownbelow. Move from right to left, noticing that A+A = A, A+B = B+A = B,and when you have B+B, you put A and "carry" the B by following thecurved line.

This is essentially the same thing that happens in ordinarybase-ten (decimal system) addition, for example when we add 9595 and405 to get 10,000.

It also illustrates how we can count arbitrarily high. All we have todo to implement "counting" is to know how to add one. And we can:given any binary number, no matter how long, we can add "B" to it byfollowing the rule described above. When needed, the rule lengthensthe number by one binary digit (bit), as when we go from 15 (BBBB) to16 (BAAAA) - just like in decimal when we go from 99 to 100. We mayor may not have invented English words for numbers withhundreds of digits in them, but they still make sense, and we know howto operate on them, just by operating on their digits.

The only drawback to binary for human (rather than computer)use is the visual monotony that stems from having only 2 symbols. Thehuman visual system can easily distinguish a much wider range ofsymbols (e.g. in English we have 26 letters, each with 2 variations:Upper and lower case). We can take advantage of this to get shorterrepresentations for our numbers. In the decimal system, we use tensymbols, conventionally written 0,1,2,3,4,5,6,7,8,9. As a result, thepositional values for the columns are powers of ten, insteadof Powersof Two. Thus for example, the number 123 in decimal means onehundred (1*10*10) plus twenty (2*10) plus three (3*1). We are so used tothis (once we have graduated from elementary school) that we think of"123" as "being" the number one hundred twenty three. Indeed ourEnglish word text parallels the decimal system very closely (with afew exceptions, like "eleven", that are common enough to get their ownword). That's why saying "10" means "two" in binary confusespeople. Hopefully now you do not mind saying that "BA" means "two" inbinary.

With the decimal system, we have more addition facts to memorize(9+6=15), but we get shorter strings of digits, which are visuallyeasier to look at. With the hexadecimal system, we carry thisidea even further: by using 16 symbols instead of 10, we can make thedigit strings even shorter, but we have yet more addition facts tomemorize. The digits 0..9 retain their ordinary meaning (when used inthe right-most column), but now we add "a" for ten, "b" for eleven,"c" for twelve, "d" for thirteen, "e" for fourteen, and "f" forfifteen. Once you get these memorized, we can do our positionalnotation trick again, this time using powers of 16. The secondcolumn (from the right) represents bundles of 16 sticks. The thirdcolumn represents groups of 16 bundles of 16 sticks, so it reallymeans 256 sticks. And so on.

For example, in hexadecimal, the symbol "a7" means "one hundred sixtyseven", or 167 in our decimal notation, because the "a" in the secondcolumn from the right tells us to take ten copies of sixteen, or 160,and the "7" in the rightmost column means to add seven more, giving167.

Because positional notation involves adding, it is notsurprising that it is easy to perform addition in any positionalsystem. That same elementary school algorithm continues to work inhexadecimal: we just need to modify the addition facts. For instance,in hexadecimal, "5+5=a", and "6+7=d", and "a+c=16", where you have tokeep in mind that in hexadecimal, "16" means 16+6, or 22. To avoidconfusion, computer programmers (who do sometimes actually need to usehexadecimal) write a "0x" in front of hexadecimal numbers todistinguish them from decimal numbers. Thus, we might write "0xa + 0xc= 0x16", and now there is no confusion: this is simply the hexadecimalway of saying that "ten plus twelve equals twenty two".

Here is the same table of numbers we used to describe the GBS system,now converted to hexadecimal.

Hexadecimal SymbolEnglish Interpretation
0x1one
0x2two
0x3three
0x4four
0x5five
0x6six
0x7seven
0xaten
0xbeleven
0xffifteen
0x14twenty
0x18twenty four
0x19twenty five
0x1atwenty six
0x1fthirty one
0x3esixty two
0x65one hundred one

As one last example, 0x123 means, in decimal: 1*256 + 2*16 + 3 = 291.Conversely, if you want to translate 291 back into hexadecimal, youneed to do a somewhat more complex process involving division: 291/256is one with a remainder of 35. Then 35/16 is two, with a remainder ofthree. So 0x123 is the answer. Isn't is nice that we have computersthat can do these kinds of mechanical translation steps for us!

I want to stress that all these multiplications and divisionsare only needed for converting between number systems. If wehad grown up learning hexadecimal instead of decimal, it would beintuitive to us, and we would not need to "convert" back to decimalin order to understand it.

Regardless of which number system we choose for writing, there is anenormously important fact about the real world lurking in thebackground: counting works. If you have five apples lined up ona table, and you count them (carefully, and without eating or droppingone), then you will always get a count of 5, regardless ofwhich end of the line you start from. If you then shuffle the order ofthe apples and recount them, you will still get a count of 5. This mayseem like a trivial observation, but it is not. There is no particularreason why the real world needs to behave in a rational, regular,systematic manner - but it does.

Said differently: as a theoretical mathematician, I can definea number system - take your pick of binary, decimal, hexadecimal, orany other - and thereby create a computational system. Ican define the operation "next", which implements counting; Ican define "addition", and all the other usual arithmeticoperations. The resulting system will be internally consistent, and Ican do all kinds of interesting things with it, BUT, none of thisnecessarily has anything to do with the real world. It is an amazingthing about the world we live in that in fact, when mathematicians dodefine new mathematical systems, they frequently turn out to haveactual predictive value in the real world. We will see an example ofthis next week when I talk about Imaginary Numbers.

At the heart of it, this is why Science works: different people, withdifferent mental models and different political agendas, cannonetheless agree that there are 5 apples on the table: not more, notless.

The real world exhibits an astonishing amount of regularityand repeatability. Not only can we count apples, but we can measuremass, position, velocity and acceleration, and summarize them via"Newtonian Mechanics", from which we can calculate how to launchrockets that can put a human on the Moon and bring them safely backhome. We can measure light, electricity and magnetism, and summarizethem via "Maxwell's Equations", from which we can calculate how tomake radios and electric lights. We can measure radioactivity and thephotoelectric effect and summarize them using "Quantum Mechanics",from which we can calculate how to make the transistors and computerchips that are central to cell phones and the Internet. In each ofthese cases, humans have been able to model the observed regularityand consistency of the real world using mathematics, then turn aroundand use that same mathematics to make predictions about thereal world, predictions that actually come true, predictions that areat the heart of all our modern technologies.


I hope you enjoyed this discussion. You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on.

As usual, please post questions, comments and other suggestions using the box below, or G-mail me directly at the address mentioned in the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can "follow" @ingThruMath on Twitter to get a 'tweet' for each new post. The Contents page has a complete list of previous articles in historical order. Now that there are starting to be a lot of articles, you may also want to use the 'Topic' and 'Search' widgets in the side-bar to find other articles of related interest. See you next time!

Genetic Diffusion Across a Geographic Barrier

To contact us Click HERE
Imagine some red flowers growing in a meadow at the base of a mountain range. Suppose one of the plants has a mutation that makes its flowers blue. If the new color makes it easier for insects to find and pollinate it, then as the years go by, the mutation will likely spread and replace the original red color in nearby plants. But will the mutation be able to spread to the far side of the mountain range? How long will it take before you find blue flowers there?

A lot of readers have enjoyed An Evolving Ecosystem Game, so today we are going to modify that agent-based simulation of an evolving ecosystem to model the diffusion of mutations across a geographic landscape with a barrier.

Let's start by describing the simulation conceptually.

We will represent land by a large grid of squares, 100 by 100.We can think of each of the 10,000 cells in our grid as having asingle plant in it. Or, for a more realism, maybe think of each cellas a square mile (or square kilometer), and assume that each cellholds a relatively homogeneous population of many plants.

In our first experiments, all the cells will hold plants. Inlater experiments, we will block off some of the cells to form abarrier, representing the inhospitable mountain range.

Each cell (plant, or "agent") has a color gene that we willrepresent as a number between zero and one, corresponding to aspectrum of colors from red to blue.

Initially, all the plants have red flowers, so we will color the cellsred. As time goes on, the plant in a given cell dies and is replacedwith a new plant grown from seed.

The color of the new plant will be inherited from its parents,which are the plant it replaced and a randomly chosenneighbor. Initially, we will assume color has no impact on fitness, sothat all neighbors are equally likely to be the parent, but later, wewill assume that colors closer to the blue end are favored bypollinating insects over colors at the red end.

Initially, I intended to let the new color be the average ofthe two parent colors, plus a small random component representingmutations that can shift it either toward red or toward blue.

However, averaging quickly wipes out the mutation, even whenit is favorable. The gene with value 1 (the initial blue plant) mixeswith one of the nearby 0's to make 0.5, which then mixes with anothernearby 0 to make 0.25, and so on. Within a few generations, everythingis back to red. This is not the situation I wanted to model.

Instead, we let the new color be inherited from one of the twoparents, chosen at random (50-50). And, with a small probability, wealso allow a mutation that shifts the gene up or down by 0.1 (thoughwe also bound it to stay between zero and one).

Implementation in Python

So - ready to try it out? As usual, today we will code up theexperiment using free, open-source software. We will once again usethe free, open-source, cross-platform programming language calledPython, which is one of the mostwidely used and easy to learn programming languages in the world.

Start by downloading and installing Python on your machine via the Python 3.2.2 Windows installer. You will also need thePyGame add-on package,usingthe PyGameWindows installer for Python 3.2, which is also free and opensource. Use the version of PyGame that most closely matches theversion of Python you are using.

My article on An Evolving Ecosystem Game has some additional pointers on getting started with Python and PyGame. I also highly recommend reading a tutorial such as the free book-length PDF file Introduction to Computer Science Using Python and PyGame.

Start Python. From the File menu, choose "New Window". This brings upa window you can use for editing Python scripts. Now copy and pastethe Python code listing below into your new editor window. Save theresulting file with a name like "flowers.py", and put it somewhere youcan find it again, such as on your desktop.

import pygameimport randomimport mathfavorable = Falsebarrier = FalsemutationRate = 0.01name = "neutral"makeMovie = Falsesize = 100 # boxes per side of gridpixels = 4 # pixels per side of boxdef breed(A):    if(A.color < 0):        return A    a = Agent(A.i, A.j)    a.color = A.color    if(random.random() < 0.5):        a.color = A.chooseNeighbor().color    if(random.random() < mutationRate):        diff = math.copysign(0.1, random.random() - 0.5)        a.color = max(0, min(1, a.color + diff))    return aclass Agent:    def chooseNeighbor(self):        a = self.randomNeighbor()        if(favorable):            b = self.randomNeighbor()            if(a.color > b.color):                return a            else:                return b        return a    def randomNeighbor(self):        i = bound(self.i + random.choice([-1,0,1]))        j = bound(self.j + random.choice([-1,0,1]))        if (i==self.i and j==self.j) or pop[i][j].color < 0:            return self.randomNeighbor() # try again        else:            return pop[i][j]    def __init__(self, i, j):        self.i = i        self.j = j        self.color = 0def bound(x):    return max(0, min(size-1, x))def rainbow(x):    z = pygame.Color(0,0,0,0)    if(x < 0):        return z    z.hsva = [x*250, 90, 90, 0]    return zslide = 1def draw(pop):    global slide    for row in pop:        for cell in row:            R = pygame.Rect((cell.i*pixels, cell.j*pixels), (pixels, pixels))            pygame.draw.rect(screen, rainbow(cell.color), R, 0)    if makeMovie:        pygame.image.save(screen, "img/slide%(n)03d.jpg" % {"n":slide})    slide += 1    pygame.display.flip()def update(oldpop):    pop = [] # new population matrix    for oldrow in oldpop:        row = []        pop.append(row)        for oldcell in oldrow:            row.append(breed(oldcell))    return pop# initialize graphics window:pygame.init()nPix = pixels*sizescreen = pygame.display.set_mode([nPix, nPix])pygame.display.set_caption("Genetic Diffusion")white = [255, 255, 255]black = [0, 0, 0]screen.fill(white)pygame.display.flip()# initialize agents:pop = [] # population matrixfor i in range(0,size):    row = []    pop.append(row)    for j in range(0,size):        row.append(Agent(i, j))if(barrier):    b = int(size/2)    for i in range(1,size-1):        pop[i][b].color = -1    for i in range(0,4):        for j in range(0,4):            pop[i][j].color = 1# main animation loop:gen = 1done = Falseclock = pygame.time.Clock()while done==False:    draw(pop)    gen += 1    clock.tick(30)    for event in pygame.event.get():        if event.type == pygame.QUIT:            done = True    pop = update(pop)pygame.image.save(screen, name+".jpg")pygame.quit()

Now, press F5 in the editor window to run the code (or go to the Runmenu if you prefer, and choose "Run Module"). The graphics windowopens and you can watch the simulation evolve. If you are using Linux,just type
python flowers.py
at the command line instead.

You should see a graphics window. Initially it will be all red, butsoon, orange and then yellow dots will appear. As time goes on, areaswith a lot of yellow will start to turn green, and if you keepwaiting for a minute or so, you will even start to see some cyan(light blue) and dark blue emerging. Ultimately, you will wind up witha picture something like this. Because of different random dice rolls,your pattern will differ in specific details, but should look prettysimilar in a qualitative sense: localized patches of all the differentcolors, all co-existing, because in this first run, weset favorable = False, so no one color is preferred.

If the code does not work for you, the first thing tocheck is whether the indentation matches what you see above:in Python, unlike in most other programming languages, indentationmatters. Python uses indentation to group statements thatlanguages like C or Java would surround with curly braces.

With this as a baseline, let's see what happens if we changethe first line of code (after the import statements) from favorable= False to favorable = True. Also, change name ="neutral" to name = "favorable". First, let's look at thecode so we know what to expect. After the "import" lines, which justload the Python packages we are using (PyGame, random numbers, and themath library), we set various global variables like the mutation rateand the size of the grid.

The breed function is central to the simulation, so we willcome back to it in a moment.

Since each agent only needs to keep track of a single number(the color gene), we don't really need the programminginfrastructure of agent based modeling here - we could just store amatrix of colors and be done with it. However, I have written the codein a more agent-based style so that it will be easier for you to use it asa starting point for making more complicated simulations.

By defining an Agent class, we create a new data type, likenumbers and strings, called "Agent". These are pretty simple agents:they have coordinates (i and j) showing where they live on the grid,and a color that starts at zero. They also have a method (function)for selecting a random neighbor, and a more interesting methodcalled chooseNeighbor. This is the function that breedwill use. It picks a random neighbor, and then if favorable istrue, it picks a second random neighbor and of the two, returns theneighbor with the largest color value.

The breed function starts with an agent called A, and checkswhether its color is negative, which means it is a black barrier cellrather than a colored plant cell. Assuming it is a real plant, webreed it by choosing a neighbor to be the second parent, theninheriting the color of one of the parents. Occasionally (based on themutation rate) we also shift the color up or down by 0.1. Finally, wereturn the new agent (plant) we have just grown from seed.

In this game, the agents are synchronized: every winter they all dieand are replaced by a whole new generation. This means we need twomatrices: one to hold the old generation and one to hold the new one,so that when we compute the parent colors for a new plant, we readfrom the old generation, not from a partial mixture of old and new.That's what the update function does, a little further down inthe code.

We also have a rainbow function for mapping the zero toone color gene into an actual Color on screen. And a drawfunction for re-drawing the screen once per new generation.

The rest of the code just puts an initial plant (agent) ineach cell of the grid, then possibly adds a barrier (mountain range),then finally runs the main loop, in which we repeatedly compute a newgeneration of plants and draw them on the screen. When you finallyclick the Close box on the graphics window, we save the final image toa JPG graphics file for future reference. For instance, we made onecalled "neutral.jpg" from the first run. When you re-run the code, youwill create "favorable.jpg", as shown below.

Actually, I stopped the simulation early - after a few moreseconds the screen would be solid dark blue. Why the dramaticdifference from the first experiment? Clearly, now that colors closeto blue are favored over colors close to red, they tend to spread. Wehave broken the symmetry, so the various colors no longer coexist inequilibrium with each other.

For out third and final experiment, set up the first few lines likethis:
favorable = Truebarrier = TruemutationRate = 0name = "barrier"
I also set makeMovie to True, and created a folder called "img"in the same folder where I have the code. This makes the code write aJPG image file for each generation, which I then combined togetherinto an animated movie. For anyone not actually running the code, youcan play the following movie to see what happens.

As you can see, we are now adding the mountain range barrier (thehorizontal black line in the middle of the picture). Up in the topleft, we create an initial mutation (in blue). In order to focus onthe spread of the initial mutation, we set the mutation rate to zeroso there will not be any others: the only way to get blue squares onthe far side of the mountain will be if the original mutation spreads.

While the screen does ultimately fill with blue, the path to gettingthere is quite different from our second experiment. This time, thegeographical barrier plays a key role in breaking the symmetry.

In case you cannot see the movie (modern browsers like Firefox shouldwork fine, but older ones from Microsoft may fail), here are threepictures from early, middle and late in the simulation, to give youthe idea.

The first time I created an animated movie from the results of asimulation was back in the early 1990's. At that time, I had to fly toCalTech inPasadena, California to use a super-computer lab with highlyspecialized (and highly expensive) equipment. The process took acouple of days (not to mention thousands of dollars). Two decadeslater, we can all do this sort of thing in seconds on our ordinaryhome computers!

I have fond memories of that trip to CalTech. One verymemorable highlight was wandering the wonderful10-acre DesertGarden at the nearby Huntington Library. If you ever have theopportunity, I highly recommend it. Coming originally from moistnorthern forests, I found it an extremely moving experience to be inthe dry desert, out of sight and sound of civilization, walking pasttowering, majestic cacti. Like seeing the Milky Way galaxy on a darknight in the deep woods, it is an experience that helps you understandyour place in the universe.

I hope you enjoyed this discussion. I encourage you to try your ownexperiments. If you change something and it fails to run, read theerror message; it often gives a good clue about what line the problemis in (usually the line you just changed), and indentation is oftenthe issue. It is best to use an editor that understands Python, suchas the one built into IDLE, since it can color code the text and helpmanage indentation for you. For instance, text following a # isignored, so we can put comments in the code to remind us whatit does.

You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on.

As usual, please post questions, comments and other suggestions using the box below, or G-mail me directly at the address mentioned in the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can follow @ingThruMath on Twitter to get a 'tweet' for each new post. The Contents page has a complete list of previous articles in historical order. You may also want to use the 'Topic' and 'Search' widgets in the side-bar to find other articles of related interest. See you next time!

Family Genetics

To contact us Click HERE
I recently read a very interesting article, The Incredible Expanding Adventures of the X Chromosomeby Christopher Badcock, Ph.Din Psychology Today.

The article explores some of the implications of the fact that fatherspass their sex chromosome undiluted to their offspring, whereasin all other situations, the offspring receives a random mix of genesfrom both parents.

The article generated considerable discussion within our household.I decided that a picture would help make things clearer. It turns outto be rather subtle: I had to re-read the article carefully and re-draw mypicture several times.

Here's how it works. Each person has 46 chromosomes arranged in 23pairs. We will ignore the 22 ordinary pairs and concentrate on thepair that determines sex. Everyone has at least one X chromosome;women have two, while men have one X and one Y. The X and Ychromosomes are quite different in length, so during reproduction theydo not inter-combine the way ordinary pairs do. In particular, themother's two X chromosomes get randomly mixed, but the father's X andY stay "undiluted": each offspring either gets the father's X,becoming a girl, or the father's Y, becoming a boy.

This apparently simple observation turns out to make for quitecomplicated mixing patterns over a couple of generations. ThePsychology Today article gives a nice text example, but lacks adiagram. Here is my attempt to fill that gap.

Let's start with your four grandparents. Assuming they are allunrelated individuals, their various X and Y chromosomes are alldifferent, indicated here by eight distinct colors. There are two Ychromosomes in the mix, from your two grandfathers: pink for thepaternal and orange for the maternal. There are also six Xchromosomes: one each for the grandfathers (purple and cyan,respectively), and two each for the grandmothers (red/blue on thepaternal side, yellow/green on the maternal side).

Now let's move to your parents. Your father, since he is male,inherited your paternal grandfather's "pink" Y. Like everyone, he getsa mix of his two maternal X chromosomes, shown here as a "red/blueflag". This mix is different for each person, so for his brothers andsisters (your paternal aunts and uncles) I have rotated the red/blueflag to make it look different. For example, your uncle is not identicalyour father, even though they share the pink Y.

You mother, since she is female, inherited your maternal grandfather's"cyan" X. Like everyone, she gets a mix of her two maternal Xchromosomes, shown here as a "yellow/green flag". Once again, this mixis different for each person, so for her brothers and sisters (yourmaternal aunts and uncles) I have rotated the yellow/green flag tomake it look different. For example, your aunt is not identicalyour mother, even though they share the cyan X.

Now it is your turn. If you are male, you will have inherited yourfather's pink Y, plus an X that is a mix of your mother's two Xchromosomes. Note that this new X is half cyan (from your maternalgrandfather), and completely lacks red, blue and purple: you havereceived contributions from only 5 of the original 8 colors. Inparticular, you have no sex chromosomes at all from your paternalgrandmother. Of course, you do share other genes of hers, on the other22 chromosome pairs - it is just the 23rd pair we are discussing here.

Similarly, if you are female, you get the red/blue X from your father,plus an X that mixes your mother's two X chromosomes. Again, I rotated thecyan/yellow/green flag to remind us that your mix is different fromthat of your siblings. Nonetheless, it still is half from yourmaternal grandfather (cyan), just as for your other siblings. Overall,you have contributions from 6 colors: you only lack the pink Y andpurple X from your paternal grandfather. Of course, you share othergenes with him, on the other 22 chromosome pairs.

That's an interesting asymmetry right there: girls have 6 of the 8colors, boys just have 5.

Also, your brothers have the same pink Y as your father andpaternal grandfather, undiluted (except for mutations) stretching backfor many generations. In contrast, although your sisters received yourfather's undiluted X chromosome, that X itself is a random mix of yourpaternal grandmother's X chromosomes, so it is not the sameas any of your grandparent's individual X chromosomes.

It's also interesting that boys share none of their colors withtheir paternal aunts. Girls, in contrast, do have some color overlap withtheir maternal uncles.

We also see that none of the grandchildren inherit any portionof the purple X belonging to the paternal grandfather. Similarly, none of the grandchildren inherit any portion of the orange Ybelonging to the maternal grandfather.

As a mathematician, I have focused here on telling a simplestory about the way the patterns change, assuming the simple rule thatfathers pass their sex chromosome undiluted to their offspring,whereas mothers pass on a random mix of theirs. In fact, an expert ingenetics would point out that real genetics is more complex than ourpicture suggests. Biology is not so "black and white" as math is!Just to hint at some of the complexities:
  • Some genes on the Y chromosome may actually mix with some on the X, so not all of them are "sex-linked" in the manner described above.
  • Some people have mutations that result in XY females, XX males, XYY males, and even other combinations.
  • In females, some cells inactivate (most of) the mother's X and other cells inactivate (most of) the father's X. This explains the phenomenon of "calico cats", where a female (XX) cat has patches of different colored hair depending on which X is expressed in that patch; as a result, there are no male calico cats.
  • The Y chromosome is much shorter than the X. That means each Xcarries lots more genes than any Y does, and men, in particular, havemany un-paired genes, where their Y lacks a counterpart to a gene onthe X.

Genetic disorders (mutations) on the X chromosome often show up inmales, because the Y chromosome does not have corresponding genes thatcould compensate for the defect. For instance, red-green colorblindness in men is actually X-linked, and hence is inherited fromtheir mothers, rather than being a defect in their father's Y.The Psychology Today article draws some interesting conclusions fromthis about intelligence and other characteristics that have somesex-linked components. You should read the original article to learn more; note that the print version in the magazine islonger and has more examples than the version available on the web.

I hope you found this interesting. You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on.

As usual, please post questions, comments and other suggestions using the box below, or G-mail me directly at the address mentioned in the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can follow @ingThruMath on Twitter to get a 'tweet' for each new post. The Contents page has a complete list of previous articles in historical order. You may also want to use the 'Topic' and 'Search' widgets in the side-bar to find other articles of related interest. See you next time!

Why are there exactly five Platonic solids?

To contact us Click HERE
You know what a cube is: a solid three-dimensional object withsix identical sides, each of which is a square. The cube is one of thefive Platonic solids: highly symmetrical shapes known at leastas far back as the ancient Greeks. In addition to the six-sided cube,we have the four-sided tetrahedron, theeight-sided octahedron, the twelve-sided dodecahedron,and the twenty-sided icosahedron. Their common feature is thatall their sides are identical. For example, all four sides of thetetrahedron are equilateral triangles, and all twelve sides of thedodecahedron are regular pentagons. Since all of their faces areidentical, they are all well suited to use as dice. Many games usecubes (common six-sided dice) since they are easy to manufacturer andship, but the popular (at least when I was growing up, before computergames appeared) fantasy role-playing gameDungeons and Dragons uses all five kinds, as shown below.

It has long been known that these five are the only shapes withthis sort of complete symmetry. Today I want to discuss the question"Why are there exactly five Platonic solids?" Why not more, or fewer?Why these particular five?

First, let's be clear about the definition of a Platonic solid. It isa geometric shape like the ones pictured above, all of whose faces areidentical regular polygons. What's a "polygon", and what is "regular"?

Well, start by drawing some dots on a piece of paper. We will refer tothem as vertices. Now connect the dots with straight lines(called edges) that do not cross, to form a closedshape. That's a polygon. For example, atriangle has 3 vertices and 3 edges.

A regular polygon is one for which all the edges are the samelength, and all the angles between the edges are the same. Forinstance, a triangle with all three edges of equal length is regular;we call it an equilateral triangle. If a four sided polygon isregular, we call it a square. If a five sided one is regular,we call it a regular pentagon.

So, to make a Platonic solid, you pick a regular polygon: either anequilateral triangle, a square, or a regular pentagon. You make abunch of copies of the polygon; we call these the faces.Finally, you glue them together along adjacent edges to form a solidthree-dimensional shape.

However, if you try this with regular hexagons (6-sided polygons), youwill find that you cannot glue them into a solid "dice-like" shape. Itjust does not work: the angles are not right. This provides a clue tofiguring out why there are only five possibilities.

Imagine we draw three squares together as shown in (A) below, then cutout the L-shape and fold on the edges between colors. The three facesmeet at a common vertex, forming half of a cube. But if we hadincluded the fourth square, as in (B), we could not achieve this: anyattempt to fold the edges down would not fit together.

Let's see if we can generalize this observation, so we can characterizewhich shapes will and will not fit together.

If we draw several identical regular polygons around a common vertex,like the three squares above, we have to leave a gap, like in (A), ifwe want to be able to fold it into a three-dimensional shape. Sincesquares have 90 degree angles between adjacent edges, four squarestake up all 360 degrees: no gap. We are limited to using threesquares, and we wind up creating the cube.

Two squares are insufficient, since we need three to form a threedimensional corner. So the cube is the only possible Platonicsolid built from squares.

Equilateral triangles have an interior angle of 60 degrees, so six ofthem completely fill the space around a vertex, just like four squaresdo. So, we are limited to using 3, 4, or 5 triangles around a vertex,if we want to build a Platonic solid. These give the tetrahedron,octahedron, and icosahedron, respectively.

With just three triangles around each vertex, the tetrahedron is very"spiky": the points are sharp, because we have to fold the papertightly to glue the edges together. Try it yourself with a piece ofpaper: print the picture above, cut out the three-triangle strip, foldwhere the colors change, and try to bend it to make a threedimensional corner. In contrast, when we make the icosahedron, with five trianglesaround each vertex we do not need to fold the paper as sharply, so thecorners are not as sharp. In fact, the icosahedron is much more like asphere (smooth all over), in comparison to the tetrahedron.

We can also use regular pentagons. The angle between adjacentedges is 3*180/5 = 108 degrees, so we can fit three together, leavinga gap, much like the three squares in (A), but four will not fit. Sothere is only one possible Platonic solid made from pentagons, namelythe dodecahedron.

Finally, if we try something with 6 or more sides, like ahexagon, the interior angles will be 120 degrees (for hexagons) ormore (for 7 or more sided shapes), so even three will fill or overfillthe 360 degrees available, with no room for a gap. So there are noPlatonic solids with faces of six or more sides. That means the five wefound are all of them: there are no others.

It turns out that there is an entirely different approach wecould have used. The analysis of angles we discussed above was knownto Euclid around 300 B.C., but this second approach is much morerecent. Various mathematicians, including Descartes in 1639, Euler in1751, and Cauchy in 1811, contributed to the discovery that for any"sphere-like polyhedron" (i.e. any geometric solid built frompolygons, without holes, such as the Platonic solids), there is arelationship between the number of vertices V, edges E, and faces F,often called Euler's formula:

V - E + F = 2

Start with a cube, for example. It has V=8 corners, E=12edges, and F=6 faces; 8-12+6 = 2 as claimed. Now imagine we take oneof the square faces of the cube and cut it into 2 triangles, asillustrated below:

What does this do? The number of vertices does not change, but thenumber of faces goes up by one, as does the number of edges. So theoverall value of V-E+F does not change, since the extra face cancelsthe extra edge. A similar argument shows that no matter how youadd or remove vertices, faces and edges, the value of V-E+F does notchange, provided you do not cut holes through the solid (i.e. turn itinto a doughnut or inner-tube shape).

Using Euler's formula, we can classify the Platonic solidswithout having to calculate with angles. Suppose the Platonic solid ismade from n-sided polygons, and suppose k of them meet at everyvertex. For instance, the cube has k=3 polygons around each vertex,and each of the polygons is a square, with n=4 edges.

If there are F faces, then before we glue them together to build thesolid, we have nF edges and nF vertices. We glue edges in pairs, sothere will be E = nF/2 edges in the final solid. We glue k verticestogether at each corner, so there will be V = nF/k overall in thefinal solid. By Euler's formula, we have

nF/k - nF/2 + F = 2.

Dividing by twice the number of edges (2E) gives

1/k + 1/n = 1/E + 1/2.

Since E>0, we know 1/k + 1/n must exceed 1/2. We also know that n isat least 3 (triangles have 3 sides, other polygons have more) and thatk is at least 3 (we need at least three faces per vertex to get athree dimensional solid).

What whole numbers n and k, each at least 3, satisfy the constraintthat 1/k + 1/n must exceed 1/2? We can list them out.

When k=3, 1/n must exceed 1/6, so n can only be 3, 4, or 5.

When k=4, 1/n must exceed 1/4, so n can only be 3.

When k = 5, 1/n must exceed 0.3, so n can only be 3.

And if k is 6 or more, 1/n must exceed 1/3, so there are no solutionsfor n (since n must be at least 3).

So once again, we have exactly five solutions, no more and noless, and once again they correspond exactly to the known Platonicsolids. The tetrahedron has k=3 triangles (n=3) meeting at eachvertex. The cube has k=3 squares (n=4) and the dodecahedron has k=3pentagons (n=5) meeting at each vertex. The octahedron has k=4triangles (n=3) at each vertex, and the icosahedron has k=5 triangles(n=3) at each vertex.

There are other interesting questions you can ask aboutgeometric shapes. For instance, someday I intend to do an article on"tessellation", which is the process of covering a sheet of paper witha repeating grid of shapes (not just squares, but triangles andhexagons work too). One can also ask more exotic questions, such as:how many four-dimensional Platonic "hypersolids" are there? But thatrequires being able to visualize geometry in 4 dimensions, which istoo tricky to delve into today.

If you enjoyed this article, you might also like Algorithmic Art. Or perhaps Counting and Number Systems.

I hope you found this interesting. You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on.As usual, please post questions, comments and other suggestions using the box below, or G-mail me directly at the address mentioned in the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can follow @ingThruMath on Twitter to get a 'tweet' for each new post. The Contents page has a complete list of previous articles in historical order. You may also want to use the 'Topic' and 'Search' widgets in the side-bar to find other articles of related interest. See you next time!

Understanding Life Expectancy

To contact us Click HERE
Demographics are destiny. As Earth reaches the 7 billion human population milestone it is more important than ever that we understand how to model the future: how fast is our population growing? What, if any, resource constraints will we face in the future as world populations continue to grow? How will the shifting balance between young and old people play out? What will be the economic, social, political and military consequences of shifting demographic patterns?

Predicting the future is hard. Nevertheless, in studying demographic trends, we have one big advantage over most social sciences: what physicists would call a conservation law. Put simply: if you are age 40 today, then in ten years, either you will be age 50, or you will be dead.

We are going to combine this fact with real data about the US from the Centers for Disease Control and Prevention (CDC) to make some graphs showing life expectancy as a function of age. Since some people will live more or fewer years than expected, we will also look at the distribution of mortality outcomes. Finally, we will also look at how life expectancy could change if medical technology improves mortality rates by 1% each year.

If we try to model the impact of advertising on shoppers, or thefuture direction of the stock market, or who will win the nextelection, we face the difficulty of predicting human behavior in theface of all our irrationality and eccentricity. But when it comes todemographics, we have an inescapable mathematical law:
population in future equals   population now plus births minus deaths

At least, that works for the total population of the Earth. If we wantto track population by country, we have to add a net migration term aswell, since not everyone stays permanently in the country of theirbirth. But we won't bother about migration today.

This fundamental law means that we know the correct mathematicalequations for modeling demographics; all that is needed is to measurethe parameters (rates of fertility and mortality). Of course, we stilldon't know the future: will people choose to have more children, orfewer, as time goes on? Will a cure for cancer result in people livinglonger, or will an epidemic kill a lot of people? But at least we canhave an informed discussion around the probability and magnitude ofthese kinds of scenarios, since the underlying math gives us theframework within which to work.

To do the calculations, we have to be a little careful aboutdefinitions. If you want, you can skip forward to the graphs, but theywill be easier to interpret if you read throughthe next few paragraphs.

We are interested in predicting the future, so we need todefine metrics that are likely to be relatively constant overtime. Suppose we track the total number of deaths in the US eachyear. This number will probably grow from year to year, even thoughpeople are living longer. Why? Simple: the US population is alsogrowing. Suppose a constant fraction of the population dies each year(e.g. 1%). As the total population P(t) grows with time, so does thetotal number of deaths D(t) = 0.01*P(t). So the ratio D(t)/P(t) is abetter metric than plain D(t), since it controls for the size of thepopulation.

However, even D/P is not that useful for forecasting,because it is not going to stay constant over time. Humanstypically live for many decades (how many varies by country), becausedeath rates are low for young people (except for infant mortality) andhigh for old people. So, if the mix of old and young changes,D/P will also change.

To illustrate, suppose the death rate is zero for people under80, and is 10% per year for people over 80. Suppose 10% of ourpopulation is over 80. Then D/P will be 0.9*0+0.1*0.1 = 0.01, or onepercent of the total population per year. Now suppose after a fewdecades, the percent of people in the population who are over 80changes to 20%. Then D/P will be 0.8*0+0.2*0.1 = 0.02, or two percentof the population per year.

Since D/P depends on the age distribution, a better model willtrack the D/P ratio separately for each age group; that way we canalso control for shifts in mix. For instance, we will keep track of D,P, and the D/P ratio for 30-year old persons, and separately also for31-year old persons, and so forth.

There is a further complication. The death rate for a given ageactually changes over time, as medical technology improves. A personborn in 1900 was 30 years old in 1930, before antibiotics were common,so their chance of dying during their 30th year was higher than thatof a person born in 1950, whose 30th year was 1980.

Demographers report death rates by age in what arecalled mortality tables. Since death rates vary as technologyimproves, we have to make clear whether the mortality rates we collectin the table refer to one cohort (e.g. the people born in 1900,followed over their lifetimes) or to one period (e.g. thedeaths that occurred during the year 2007, which impacted people froma whole spectrum of birth years).

Constructing a cohort table requires taking data over ahundred year period, and you have to wait until all of the members aredead before you are finished. That makes period tables much morepractical.

Notice, though, that period tables come with the implicitassumption that their rates stay constant from here on out. Thisassumption can be unrealistic. Imagine a baby born in 2007. Theyexperience the infant mortality conditions of the table accuratelyenough, but as they age, they ought to experience less mortality thanthe table suggests. That's because when they are 30, for example, theyear will be 2037, and (hopefully) improved medical care will reducethe death rates for 30 year old people compared with the value in thetable, which reflects the experience of people who were 30 in2007. Similarly, when this 2007 baby is 80, the year will be 2087, not2007, and so the mortality rate ought to reflect another 80 years ofprogress. Thus, the period table shows the life expectancy that thebaby of 2007 would experience if medical care never improvesbeyond the 2007 level.

The bottom line is that the period life table gives us apractical way to compare death rates by age at one moment in time,even though it cannot accurately predict what will happen to yourcohort (let alone you personally), in some distant futureyear. Comparing period tables constructed in different calendar years(or for different countries or races or genders) also gives us astandardized way to compare and track progress.

Demographers measure mortality rates (and also fertilityrates) subdivided by country, age, gender, and sometimes additionalfactors like race. They also have to track immigration, because peoplecan change countries over time. For today though, we are going tolook just at mortality, just in the United States, and withoutsubdividing by race or gender. We will leave fertility and otherelements of a complete model to a future article.

The Centers forDisease Control and Prevention (CDC) publishes "Vital Statistics" eachyear, including aso-called LEWK3 table; we will work with the most recent one, from 2007.This is a period table. From it, we can calculate all sorts ofinteresting metrics, such as life expectancy.

Here are the first few rows in the table - the spreadsheet, which youcan download from the CDC, actually goes all the way to age 100.

Let's be sure we are clear about what each column means. From the CDC documentation file nvsr58_10.pdf page 2:
  • The row labeled Age 4-5 refers to what happens to people during the year between their 4th and 5th birthdays, i.e. while they are 4 years old.
  • The probability of dying column (q) is the probability of the person dying during that year of age. Observe that it is much higher (more than 10 times) for age 0-1 than for subsequent years: this is why infant mortality is such a big concern.
The documentation is a bit ambiguous: is q[A] the probability of deathat age A for all people in the original cohort, or just for people whohave attained age A? It must be the latter, since column (q) does notsum to one, as it would need to if these probabilities applied to thewhole cohort. But since all the other columns in the spreadsheet arecalculated from column (q), we can verify this interpretation byreplicating the calculations ourselves.

The number surviving column (l) starts at 100,000 people and goes downover time as the people age. Here is where we blur the distinctionbetween period and cohort. The probabilities (q) are period deathrates, e.g. in 2007 in the US, 0.0460% of children age 1-2 died, and0.0176% of children age 4-5 died. To calculate column (l),we pretend that these rates will not change in thefuture. Hence, our hypothetical cohort of 100,000 people willexperience the infant mortality of the 2007 cohort, the age 1-2 mortalityof the 2006 cohort, the age 2-3 mortality of the 2005 cohort, and soon, since all of those rates were observed in 2007 and hence recordedin our table. As noted above, if medical technology continues toimprove, then for the real 2007 cohort, by the time they reachany given age, their actual mortality should be lower than given inthe table, so more will survive longer than predicted in the table.

But assuming medical care no longer improves, we can followour hypothetical cohort forward and observe what happens. Multiplying100,000 by 0.006761 gives 676 infant deaths, leaving99,324. Multiplying 99,324 by 0.000460 gives 46 deaths during age 1-2,leaving 99,278 survivors. And so forth. If you download thespreadsheet, you can work with exact numbers and avoid round-offerror. Following this pattern you can verify the (l) and(d) numbers from (q).

Since there are one hundred independent numbers in the table(the contents of the (q) column), the best way to understand the tableis via a graph. However, people often prefer having a single numberthat summarizes a lot of the information in the table. There arevarious numbers you could compute for this purpose, such as theaverage death rate, but the one most commonly reported iscalled life expectancy at birth. This is the number 77.9 foundin the age 0-1 row of column (e).

Life expectancy means the expected, or average, number ofadditional years a person in our hypothetical cohort will live. How shall we calculate it? Among people who reached age x, some die inthe next year, during their x-to-(x+1) year. Some die the followingyear. Some keep living a long time further. We have to weight each ofthese outcomes by the corresponding probability, and then add them upto find the expected value.

The two columns labeled (L) and (T) help with this calculation.Except for the first and last rows in the table, the (L) column showsthe average population during that year, assuming people who die do soat a uniform rate throughout the year. For example, the (L) value forage 3-4 is 99,239, which is half way between the surviving population(l) at the start of the year (99,250) and at the end of the year(99,228). The first and last values have been modified based oninformation not shown to us, namely the age in months at which infantsdied, or the age beyond 100 that people survived to. Apparently theyassumed people live exactly two years past 100.At any rate, (L) represents the total number of person-yearslived during this age range: one year for everyone who survives theyear, plus on average half a year for those who die during the year.

The (T) column is calculated from the bottom end of the tablebackward to the top, by summing the (L) column. For instance, thevalue of (T) in the 0-1 age row is the sum of the entire (L) column.Since (L) represents people-years lived during one age row, (T)represents the total people-years lived during and after one age row.Since those total people-years belong to the survivors listed incolumn (l), if we divide (T) by (l) we get the life expectancy beyondthe current age.

Personally, I find it easier to think about this definitionthe other way around: take the probability-weighted sum of the extranumber of years people live. The tricky part is that the probabilitiesare not given by the (q) column, but rather, by the (d) column fromage x onward, divided by the (l) column at age x. To see this, notethat the people in cell l[x] die according to the values in d[x],d[x+1], and so forth, which therefore sum to l[x]. So, for example,starting at age x=10 (meaning row 10-11), if we multiply column (d) by(age-9.5), sum, and divide by l[x], we get 68.6, as expected from thespreadsheet. The spreadsheet formulais =SUMPRODUCT(H15:H105,D15:D105)/C15 after firstputting =0.5+ROW()-15 in column H.A little algebra should convince you the two approaches give the sameanswer.

So far, all we have done is verify the definition of lifeexpectancy. Now for the much more interesting part. Let's draw somegraphs, and then let's ask how these concepts may change as time goesforward. In other words, if medical technology does continue to improve,how much longer might we expect to live? And, rather than just look ataverages, let's look at the distribution of outcomes. There isso much more information here than can be summarized in a singlenumber!

If you like, you can do the calculations and make the graphs yourselfusing a spreadsheet. However, I will show how to do it using the freehigh-quality open-sourcestatistical programming language R. You can follow along bydownloading your own completely free copy of Rfrom The Comprehensive R Archive Network. Look in my earlier post onKoch Snowflakes for some background on why R is a good choice, or look in System Dynamics: Feedback Models for more on population modeling using R, or look in How Do We Know? for an even simpler introduction to R.

The first thing we need to do is use our spreadsheet program toconvert the CDC table, which was designed for human readers, into aplain tab-delimited text file, with numbers in the age column, like this:

Save this in afile called "us-period2007-life-table.tab". Now we can close off thespreadsheet and do the rest with R.

Here's the R code we will use. The first line reads the tableand turns it into a "dataframe", which allows us to refer to the columns by name. The secondline prints a summary of each column so we can check that it read thefile correctly. Next, we define a function for calculating the (l),(d) and (e) columns from the (q) column; this will let us experimentwith modifying the mortality rates. The rest of the code makes thegraphs that we will discuss below.

mort <- read.table('us-period2007-life-table.tab',  header=TRUE, sep='\t')print(summary(mort))N <- length(mort$q)lifeExp <- function(d) {  y <- 1:N  p <- 0  np <- d[N]  for(i in N:1) {    p <- p + d[i]    np <- np + (i-0.5)*d[i]    y[i] <- np/p-i+1  }  y}calc <- function(q) {  l <- 0*q + 100000 # initial cohort size  d <- 0*q  e <- 0*q  for(i in 1:N) {    d[i] <- q[i]*l[i]    l[i+1] <- l[i] - d[i]  }  list(l=l, d=d, e=lifeExp(d))}m <- calc(mort$q)## check we can reproduce 'l', 'd', and 'e' from 'q'if(m$l[N+1] != 0 ||   max(abs(m$l[1:N]-mort$l))>0.1 ||   max(abs(m$d-mort$d))>0.01 ||   max(abs(m$e-mort$e))>0.8)  stop('does not match data')## now make some graphics helper functions:graph <- function(name, ylabel, y, col='black', extra=0) {  png(paste(name,'.png',sep=''), 800, 500)  par(mar=c(5, 5, 1, 1), cex=1.5, lwd=2)  n <- length(y)  yy <- y[1:(n-1)] # drop last point  plot(c(0,n-2), c(0,max(yy)+extra),  type='n', xlab='Age', ylab=ylabel)  add(y, col)}add <- function(y, col='black') {  n <- length(y)  lines(0:(n-2), y[1:(n-1)], col=col)}## now make the plots:graph('le', 'Life Expectancy', lifeExp(mort$d))dev.off()graph('mr', 'Mortality Rate (%/year)',        100*mort$q[1:(N-1)])dev.off()graph('sv', 'Surviors (% of cohort)', mort$l/1e3)dev.off()## modify mortality rates up or down 50%:m <- calc(mort$q * 0.5)m2 <- calc(mort$q * 1.5)graph('le1', 'Life Expectancy', lifeExp(m$d), 'blue')add(lifeExp(mort$d))add(lifeExp(m2$d), 'red')dev.off()graph('sv1', 'Surviors (% of cohort)', m$l/1e3, 'blue')add(mort$l/1e3)add(m2$l/1e3, 'red')dev.off()## make a cohort graph following 2007 or 1957 people,## assuming a 1%/year improvement after 2007m <- calc(mort$q * (0.99^(1:N)))m2 <- calc(mort$q * (0.99^((1:N)-50)))graph('le2', 'Life Expectancy', lifeExp(m$d), 'green')add(lifeExp(mort$d))L <- lifeExp(m2$d)lines(50:99, L[50:99], col='blue')dev.off()## age of death probabilities for the 1957 cohort, ## again assuming a 1%/year improvement after 2007graph('death', 'Probability of Death by Age',      m$d/m$l[1], col="black", 0.003)lines(50:99, m2$d[50:99]/m2$l[50], col="blue")dev.off()

First, the life expectancy graph, column (e) in the dataset. This shows that at birth, life expectancy is about 78, falling roughly linearly with age until around 60, after which it starts to tail off. Since not everyone dies at age 78, the graph has to tail off: it cannot go zero, let alone be negative. For example, at 80, life expectancy is about 9 more years. This means that having survived all the way to 80, you will, on average, live 9 more years - even though 80 is already past the original life expectancy at birth of 78. Of course, "you" means a member of the hypothetical cohort born in 2007 experiencing no further improvement in medical care beyond 2007 levels, not the real "you". If you are already 80, today when you read this, then these numbers are representative for your generation, but if you are 20, you can hope that 60 years of progress will make the numbers for your generation better, assuming you actually make it to age 80 yourself.

Next, we plot the individual mortality rates for each age, column (q)in the dataset. The rate for age 100 is 100%, but that is an artifactof ending the table there, so I have suppressed it to enlarge thevertical scale. You can see the blip at zero for "infant mortality",followed by almost zero death rates until people reach age 60 or so.

Finally, we look at the percent of the cohort surviving to a givenage, column (l)in the dataset. Aside from the infant mortality blip at the start,this decays very slowly until around age 60, after which itaccelerates for a while, and then tails off, since like lifeexpectancy, it can never go negative.

Visually, it looks like the median age of death (the age where the survivalcurve is at 50% of the population) is around 80, which makes sense:the median will be fairly close to the mean (life expectancy), thoughnot identical.

Why did we bother figuring out how to calculate these columns? Why notjust graph them directly in our spreadsheet? Well, now that we knowhow to do it, we can change things. In particular, we can experimentwith reducing the death rates to see how much of a difference medicalimprovements might have on life expectancy.

Since we do not know the future, we will have to presentscenarios. The next two plots show life expectancy and survival curvesunder the following scenarios:
  • black: the base case shown above
  • blue: mortality rates are cut in half across the board, meaning we multiply column (q) by 0.5
  • red: mortality rates rise by 50% across the board, meaning we multiply column (q) by 1.5
Can you guess what will happen?

A fifty-percent change in mortality rates seems like it should havepronounced impact. Intuitively, we expect something like a fifty-percent change inlife expectancy! In fact, life expectancy at birth changesby only 5 years up or down, with progressively less impact for olderpeople. What is going on?

Think of it this way: until age 60, the mortality rate isalmost zero, so whether we multiply it by 0.5 or by 1.5, it is stillalmost zero. Only when people are 70 or 80 or 90 are the rates highenough that 50% up or down is a big change. As a result, lifeexpectancy at birth still reaches into the 70's, even for the redline, and does not get far into the 80's, even for the blue line. Thatmeans that changing life expectancy by even one year is hard. Thatsuggests that the differences in life expectancy betweenindustrialized countries and developing countries will be difficult toreduce without massive improvements in health care;see GapminderWorld Map (2010) for a well-drawn diagram showing the correlationbetween life expectancy and economic progress.

Similarly, if we look at the survival graph, most people live to their70's or 80's under all three scenarios. For ages below 60, thedifferences are small. "Almost zero" mortality rates do add up overtime, but it takes half a century or so to see it.

So far, these graphs reflect the period life table - they showwhat would happen to the cohort of babies born in 2007 if medicalprogress freezes at 2007 levels. Let's try to make some cohortgraphs, predicting what will actually happen for those babies, as wellas for the cohort born in 1957, which is always 50 years older thanthe 2007 cohort. To draw these graphs, we have to make an assumptionabout how fast mortality rates will improve in the future. For simplicity- since I have no real data on this - let us assume that the (q)values will fall by one percent (multiply by 0.99) with each passingyear. In the next graph:
  • The black line is the usual 2007 cohort life expectancy as in the previous graphs,
  • The green line is the "actual" average number of remaining years of life for people from the 2007 cohort who have attained a given age, based on our 1%/year improvement assumption, and
  • The blue line is the "actual" for people from the 1957 cohort.
Can you guess what it will look like?

To make the green line, we multiply the mortality rate for age A by0.99^A, reflecting A years of progress since 2007.

However, to make the blue line, we multiply the mortality rate for age A by0.99^(A-50), since the blue cohort was already 50 years old in 2007.

The blue line begins at age 50, since the 1957 cohort has alreadyreached 50 as of 2007 - all these graphs are drawn from theperspective of 2007, not 2011, since there is apparently a 3 year lagin publishing the table.

We see that life expectancy at birth for the 2007 babies is really 83,up 5 from 78 years, provided medical care improves as we haveassumed. For those who survive to age 50, their expected additionalyears of life after 50 also rises by 5 years, from 32 to 37 (i.e. theycan expect to live to 87, rather than 82).

For people born in 1957, though, their expected additionalyears of life as of age 50 only increases by 2 years, from 32 to 34,i.e. they can expect to live to 84, rather than 82. That's becausetheir old age will happen relatively soon - in just 30 years - sothere will not have been time to improve medical care as much as forthe 2007 cohort, for whom old age (i.e. years with large mortalityrates) is still 70 years away.

All the graphs so far looked at averages. For any individual person, though,it is also interesting to know the distribution of additionalyears of life. After all, some people die young, while others livepast 100. Not everyone lives exactly to the average. What does the "bell curve" look like?

Turns out there is an easy answer. We just need to look at column(d), which holds the number of people from the original cohort thatdied at each year. Starting at a particular age A, we divide by thetotal number of people that reached A, namely l[A], to get theprobability distribution.

In the final chart, the black line shows the probability of death at agiven age for the 2007 cohort, and the blue line for the 1957cohort. In both cases, the calculations are as of 2007, and assume a1% improvement in mortality rates per year throughout the 21st century.

This picture is quite interesting. First of all, on the left side, theinfant mortality piece looks much more noticeable than it did backwhen we plotted the mortality rates. That's because those ratesleave out the size of the population they apply to. The population islargest at birth, and declines over time, so the high rates at agesnear 100 do not actually translate into very many people, since somany died along the way. This graph makes it more clear that infantmortality is still a very big problem even in the US.

Looking at the right-hand side, we see that the peak for the 2007cohort comes at a later age than for the 1957 cohort - again becausethe 2007 cohort has an extra 50 years of medical improvements to helpthem over the 1957 cohort. Wecalculated earlier that the 1957 cohort could expect to live to 84 onaverage; now we see that there is considerable spread around thatnumber, with significant numbers of people dying at every age between50 and 100.

In fact, we see that cutting off the life table at age 100 is a bitpremature: since it is all computerized nowadays, there is no reasonnot to extend it out to 110 or even 120 in order to provide moreinsight into how medical improvements affect the oldest people.

I hope this example encourages you to experiment: in whatever countryyou live, download the latest life tables and see what the forecastsare for you. You could also download several life tables fromdifferent years and compare them in order to estimate just how muchprogress improved medical care is making each year. Of course, pasttrends need not continue in the future, but they are at least astarting point for discussion.

I hope you found this interesting. You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on.As usual, please post questions, comments and other suggestions using the box below, or G-mail me directly at the address mentioned in the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can follow @ingThruMath on Twitter to get a 'tweet' for each new post. The Contents page has a complete list of previous articles in historical order. You may also want to use the 'Topic' and 'Search' widgets in the side-bar to find other articles of related interest. See you next time!

7 Temmuz 2012 Cumartesi

Heat Waves and Global Warming

To contact us Click HERE
A new scientific report says that extreme heat waves will likely become much more common in the future, due to global warming. For instance, the sort of extremely hot day that used to occur only once in twenty years will likely become an every-other-year event by the end of the 21st century.

This may seem rather surprising. We have heard that Earth's average temperature will likely rise over the next century, but the amount sounds modest: only a few degrees, small compared to the 30 degree (Fahrenheit) swing in temperatures between day and night, or the 60 degree swing between winter and summer.

While rigorous climate modeling is very technical, it turns out that we can get a basic understanding of the reason for the dramatic increase in extremes using some very simple calculations, which are the focus of today's post.

The Associated Press ran a story summarizing thenew report from the Intergovernmental Panel on Climate Change.The report uses careful language to clarify the strength of the evidencefor each of its conclusions. The relevant quote for us is:
It is very likely that the length, frequency and/orintensity of warm spells, or heat waves, will increase over most landareas. Based on the A1B and A2 emissions scenarios, a 1-in-20 yearhottest day is likely to become a 1-in-2 year event by the end of the21st century in most regions, except in the high latitudes of theNorthern Hemisphere, where it is likely to become a 1-in-5 yearevent.

The report covers other weather phenomena besides heat waves, but Iwant to stick to heat extremes so that we do not get sidetracked bydiscussions about the physical mechanisms that might connect highertemperatures to more precipitation, more drought, or more storms.

Here's the basic idea for understanding how extremes get magnified:
  • We start with a "bell shaped curve" representing the variety of temperatures we see in today's climate.
  • We shift the curve to the right by a small amount, representing the rise in Earth's average temperature over the next century.
  • We examine what happens to the high temperature part of the curve, i.e. to the frequency of heat extremes.
First, let's look at a picture.

I have labeled the vertical scale "hypothetical" frequency toemphasize that this is not a plot of real temperature data.It is just a "bell shaped curve" designed to illustrate the ideas.

The black line represents the original 20th century climatefor some hypothetical location, with an average temperature of 40degrees. The taller parts of the graph correspond to the more commontemperatures: in this location, most temperatures are between 0 and80, with temperatures as low as -40 or as high as +120 extremely rare.

The red line represents the climate at some future point whenthe average temperature has risen 5 degrees, to 45. In this example,we have not changed the shape of the curve at all: we have simplyshifted it five units to the right (real climate modeling alsoattempts to understand how the shape might change). This makes thetypical temperature range 5 to 85. Temperatures over 100 arerare, since both curves are close to zero near the righthand side of the chart.

Now let's now zoom in on the right-most part of the curve: the"right tail", representing the very hottest days.

Look at the blue line. Days over 110 are quite rare (about one in 12years for the black line), but they happen more than twice asoften under the red curve (about one in 5 years). So even a small shift inthe curve (5 degrees out of a typical spread of 80 is about a 6%shift) can produce a large change in extreme conditions (2.5 times asoften, in this example).

My example curves were drawn using a "normal" distribution,which is the widely used "bell-curve" shape of classical statistics,but there is no particular reason that real temperature data needs tobe shaped that way. Let's now look at some actual temperature data,and see if we can illustrate a 1-in-20 year event becoming a 1-in-2year event.

You can get real temperature readings from many sources. Somecountries make it particularly easy to download large amounts data for free.Australiahas a particularly user friendly site. I downloaded the daily maximumtemperature (degrees Celsius) in Melbourne, station number 086071, forthe years 1868 through 2010, inclusive, which provide 52,230 qualitychecked daily values with no missing values.

We will use the following R code to make the graphs. Asdiscussed in many previous posts, you can follow along by downloadingyour own completely free copy of the open source, high quality,statistical programming language R fromThe Comprehensive R Archive Network.

ds <- read.table("melbourne.tab", header=TRUE, sep="\t")dc <- 3.5 # shift in degrees Cpng('melbourne.png', 800, 500)par(mar=c(5, 5, 1, 1), cex=1.5, lwd=2)plot(density(ds$maxDayTempDegrC+dc), col="red",     xlab='Temperature (F)', ylab='Actual Frequency (%)')lines(density(ds$maxDayTempDegrC))dev.off()z <- rev(sort(ds$maxDayTempDegrC))[7] # 44.1 Cpre <- sum(ds$maxDayTempDegrC >= z)post <- sum((ds$maxDayTempDegrC+dc) >= z)print(c(pre, post, post/pre))png('melbourne2.png', 800, 500)par(mar=c(5, 5, 1, 1), cex=1.5, lwd=2)plot(density(ds$maxDayTempDegrC+dc), col="red",     xlim=c(40,50), ylim=c(0,0.006),     xlab='Temperature (F)', ylab='Actual Frequency (%)')lines(density(ds$maxDayTempDegrC))lines(c(z, z), c(0,0.005), col='blue')dev.off()

First, we plot the temperature distribution from 1868 through 2010 inblack, then overlay a shifted version of the curve in red, based onadding 3.5 degrees C to all the temperatures.

Next, we zoom in on the section above 40 Celsius (104 Fahrenheit):

The sample sizes for extreme events are by definition low,so by nature this sort of prediction isvery sensitive to the exact nature of the shift in the curve(i.e. does it also change its shape, not just slide over), as well asto where we draw the line to denote "extremes".

If we choose 42 C (108 F) as out threshold, days that hot orhotter would be expected to occur ten times more often in future thanin the past, because while only 40 past observations were that hot,400 future ones are.

If we want to look for a 1-in-20 year event in the historicaldata, we need a temperature threshold that is reached or exceeded only7 times among the 52,230 observations. That is 44.1 C in this data set(the vertical blue line in the last chart). Temperatures reach or exceed thisin the hypothetical future data set 109 times, which is 16 times moreoften, making it an almost every year occurrence!

Climate change is a complex subject, and careful data analysisand modeling require great technical sophistication. Nonetheless,these examples illustrate the simple point that while we may not knowthe exact amount of the impact, a relatively small shift in theaverage value of the temperature distribution can indeed produce verylarge increases in the relative frequency of rare events like heat waves.

I hope you enjoyed this discussion. You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on.

If you enjoyed this topic, you may also like Tracking Hurricane Irene.

Please post questions, comments and other suggestions using the box below, or G-mail me directly at the address mentioned in the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can follow @ingThruMath on Twitter, where I will tweet about each new post to this blog. The Contents page has a complete list of previous articles in historical order. Now that there are starting to be a lot of articles, you may also want to use the 'Topic' and 'Search' widgets in the side-bar to find other articles of related interest. See you next time!