# St Petersburg Paradox

One of my favorite courses in college was “Experimental Economics”, where we examine how real live humans differ from homo economicus- the perfectly logical, profit-maximizing agent of classical economics.

A good example is the ultimatum game, a 2 player game where one person has $10, and suggests a way to split that$10.  The second player can accept, and then they go on their merry way, or decline and both players get nothing.  Homo economicus would accept $0.01, since she would be better off than when she started, but most of us would be pretty insulted by such an offer. Conversely, most of us probably would not offer only$0.01, even if it is the “best play”.

Another game I am intrigued by is called the St. Petersburg paradox, apparently first posed and studied by some combination of the Bernoulli family.  The idea is that we will play a game where I flip a coin until it comes up tails.  If it comes up heads n times in a row, you win $2^n$ dollars.  The question is: how much would you pay to play this game?

A reasonable way to approach answering this is to compute that you have a 1/2 chance of winning $1, 1/4 chance of winning$2, 1/8 chance of winning $4, and so on. Hence, your “expected winnings” are $\frac{1}{2}1 + \frac{1}{4}2 + \frac{1}{8}4 + \cdots = \frac{1}{2} + \frac{1}{2} + \frac{1}{2} + \cdots$, which is a divergent sum. This means that if I offer to let you play this game for$100, a “supremely logical” person would accept.  Or for $1,000,000. Or any number. On the other hand, I know very few people who would pay even$10 to play the game. The wikipedia article above has plenty of good explanations for why, but this strikes me as a great example of where human behavior diverges from what is “optimum”.

Since I also like creating some content every once in a while, I also have a python script that plays this game using a random integer generator.  Playing the game 10,000,000 times, the average payout was $11.66, and the largest payout was$4,194,304, meaning 22 heads in a row.  This is somewhat meaningless, as we already calculated that the average payout will grow without bound as I perform more trials, but gives a feel for how much the game might be worth intuitively.

from random import randint
def game():
coin = randint(0,1)
count = 0
tot = 1
while coin != 0:
tot+=2**(count)
count +=1
coin = randint(0,1)
def lotsogame(N):
tot = 0
top = 0
for j in range(N):
n = game()
tot += n
if n>top:
top = n
return float(tot)/float(N), top
print lotsogame(1000000)

# Late to the recursion game.

Ok, confession time: every once in a while, you find out something that apparently everyone except for you knows about (I’m looking at you, pronunciation-of-rendezvous).  This time, I totally missed out on recursion in the course of teaching myself how to program.  For those in the same boat as me, this means you can call a function inside the same function.  So long as you specify a base case (and make sure you reach this base case), you can get a reasonable answer.

In the course of rectifying this, I’ve written two toy programs.  The first is very un-mathy, but essentially solves the game “TextTwist“. First I use recursion to find all possible substrings of a collection of letters, then I import an English dictionary (thank you, Python!), and check these substrings against the dictionary to come up with a list of English-language words.  See below for the print-out.

The dictionary I am using is a little bit generous with what it allows as a word, but I'd rather err on the side of being thorough in this case.

The second program I wrote was in MATLAB, and the code is a little cleaner, so I can in fact post it here.  It uses recursion to again generate the Sierpinski gasket from the flurry of fractal posts I had earlier, though this is a deterministic approach.  What it does is takes a column of points determining the vertices of a triangle (or any other shape), then for each point, returns a column determining the points of a triangle whose vertices are halfway to the previous vertex, and halfway to the next vertex.  The code may be easier to read, noticing especially that I call gasket on the third to last line:

 function [A,B] = gasket(X,Y,N) %returns the coordinates of the vertices for a Sierpinski gasket with N %iterations m = size(X,1)-1; A = []; B = []; for k = 1:m-1 A = [A,[X(cmod(k,m),:);(X(cmod(k+1,m),:)+X(cmod(k,m),:))/2;(X(cmod(k-1,m),:)+X(cmod(k,m),:))/2;X(cmod(k,m),:)]]; B = [B,[Y(cmod(k,m),:);(Y(cmod(k+1,m),:)+Y(cmod(k,m),:))/2;(Y(cmod(k-1,m),:)+Y(cmod(k,m),:))/2;Y(cmod(k,m),:)]]; end if N~=1 [A,B] = gasket(A,B,N-1); end end

The function “cmod” just changes the “mod” function so that “cmod(m,m) = m” instead of the usual behavior where “mod(m,m) = 1″.  Running this code with N = 8, we get the following image:

Noticeably cleaner than at least some of the randomly generated fractals from earlier, and at the very least faster.  After 8 iterations, there are 3^8 = 6561 triangles, and for each we store 4 pairs of numbers to generate a triangle (I am calling the “line” command, which means I need to specify that the start and end points are the same), so the total number of numbers I store is 52,488, compared to 100,000 for this much lower resolution version from earlier.  Also, I can start with any set of points I want with this program, so here is a belated Valentine for everyone:

It is somewhat hard to find interesting, mostly convex parametric curves, ok?

# Patterns in numbers

I am partial to some examples of coincidence (or not). Hint: see the left column of letters.

I tend to shy away from even faint whiffs of mysticism, numerology, and in general attaching undue weight to circumstance.  As such, I was initially a little upset when sent this link from my friend over at the fantastic Austin music blog ovrld.com.  For those who don’t wish to click, the link points out that the fraction 1/998,001 prints out the digits from 000-999, in order, except for 998.  I mean, if we look at the first 1,000,000 fractions, we’ll to find 1,000,000 decimals that might hide some patterns significant to someone – and that’s only fractions of the form 1/n.

Also, recall that if you want a repeating decimal with some pattern, there is a nice algorithm to find the associated fraction: to find the fraction equal to x = 0.19851985… (being both my day of birth and my year of birth), we notice that

10000x-x = 1985.1985… – 0.19851985… = 1985.

Then solving for leaves us with x = 1985/9999.  Here’s the thing though: this is not very pretty- the numerator is pretty big.

So I guess some kudos should go to whoever noticed this fact about 998001.  Let’s just not attach any special significance to it.

Today's filler photos come from Shiner, TX. Czech it out!

As an aside, I assume this was done by noticing that

1/98 = 0. 01 02 04 08 16 32 65 30 61 22 44 89 79 59 18 36 73 46 93 87 75 5…,

which is almost the sequence of square numbers powers of two (once you get to 128, the “1″ gets added to the “4″ in “64″, and so on), which is somewhat interesting.  Looking at higher denominators like this, you might get to

1/9801 = 0. 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 99…

Still don't know what Lady Liberty was doing in Texas...

Which might somehow suggest getting to 998001… eventually.  Finally, it turns out that this fraction is not so special, in the following way: it does not start repeating for another few thousand digits.  In particular, it goes through 2997 digits, then cycles back to the “000001002003…” pattern.

I wrote a very short Python script to play around with this sort of thing- it is a pretty quick way to perform division arbitrarily accurately, and just uses the Python function  divmod.

def long_div(p,q,digs):
#finds the first digs digits of p/q, returned as a string.  p and q must be integers.
a,r = divmod(p,q)
n = str(a)+'.'
for j in range(digs-len(str(a))):
r *=10
a,r = divmod(r,q)
n += str(a)
return n

# A function for displaying bits

A very busy week, so today’s post is just a code snippet which converts a decimal number to bits.  Notice that you can go inside the code and adjust how many bits are given to exponents and to the mantissa.  I’ve adjusted a bit to try to help the formatting, but I want to prioritize copying and pasting.  Let me know any problems/criticisms of the code!

from math import log,floor

def bitrep(mbits,ebits):
#this function writes a float N = sign*M*2**(E-offset), where sign = 0 or
#1, M ranges from 2**0 up to 2**(mbits+1)-1 and E ranges from 2**0 up to
#2**(ebits+1)-1.  The numbers sign,M,N are then converted to binary, and
#sign,M,N are then returned.  To get the bit representation of a double,
#set mbits = 52, ebits = 11.
mbits +=1
offset = 2**(ebits-1)-1
N = input('What decimal do you want the bit representation of?\n> ')

#find the sign of the number
if N >= 0:
S = '0'
else:
S = '1'

N = abs(N)
#find the size of the number's exponent
E= int(floor(log(N,2)))+offset-mbits

#convert N into the fraction n/d, where d is a power of 10.
n = str(N)
try:
d = 10**(len(n)-n.index('.')-1)
n = int(d*N)
except(ValueError):
try:
d = 10**(len(n)-n.index('e')-1)
n = int(d*N)
except(ValueError):
d = 1
n = N

#finds the mantissa of N
if offset>=E:
if d == 1:
M = N*2**(-E+offset)
else:
M,r = divmod(n*2**(-E+offset),d)
if r>=d/2:
M+=1
else:
M,r = divmod(N,2**(E-offset))
if r >= (E-offset)/2:
M+=1

#converts both the exponent and mantissa to binary
Ebin = bin(E-1)
Ebin = Ebin[Ebin.index('b')+1:]
Mbin = bin(M)[3:]

#pads each number up to the total number of bits
Ebin = Ebin.zfill(ebits)
Mbin = Mbin.zfill(mbits)
#cuts the mantissa back down to mbits of length, if necessary
Mbin = Mbin[:mbits-1]

#print out with the mantissa and exponent, then carrying out the
#calculation with python, then in bits
print N,'\t"="',M,'* (2**[',E,'-',offset,'])\n','\t = ',M*(2**(E-offset)),'\n',S,Ebin,Mbin

bitrep(52,11)

# Getting started with Python, on a mac

I think the default Mac terminal colors are different from this, but you can change them in "preferences"

This is intended as a quick guide, since one of the things that got me into programming was realizing how quickly I could start on my Apple laptop without installing anything.  As a practice first project, we’ll first compute the first 10 Fibonacci numbers, then write a function that will do this.
1. Finding Python.  The easiest way to get started is to just hit “cmd+space” (this opens spotlight: my favorite Mac feature), then type “terminal”, and open the terminal.  After that, it is as easy as typing python, and you have either a powerful calculator, or a way to start writing scripts.
2. Fibonacci. As a nice first calculation, we’ll save two variables to initialize the Fibonacci sequence, then run a loop 99 times (I had to play around a bit to see how the loop count corresponded to where in the Fibonacci sequence I was). See the image for this section.  Note that the spaces are very important.  You need to end lines with a colon, and start lines after the colon by indenting.  See below for more proper indenting.

Wikipedia claims the first Fibonacci number is 0, which I am ignoring for this post.

3. Function. Now to write a helper function that I can use to calculate many such numbers, I will open a text file (just in TextEdit), type in the code below, and save it as Fib.py:

 def Fib(n): #calculate the nth digit of the Fibonacci sequence a,b = 1,0 for j in range(n): a,b = a+b,a return a

for j in range(10): print Fib(j)

You can also see this post in progress behind the semi-opaque terminal!

4. Running a script. Now go back to the terminal window, hit “ctrl+d” to get out of Python, and navigate to the folder where you saved Fib.py. Mine was in “colcarroll/mystuff/blog”. Once there, just type “python Fib.py“, and the terminal will print out anything you asked to be printed out. Above, I print out the first 10 Fibonacci numbers.

Cheers! Now go find a proper introduction. Or use Project Euler.

# Floating point, continued

As a followup to yesterday’s post on a floating point problem I encountered using the arctangent, I want to discuss how floating point numbers are stored in computers.  Apologies in advance to my misrepresenting a lot of this (and corrections are, as always, welcome), but it is helpful for me to work it out.  If you want to just get straight to experimenting, I am writing a Python script which will be posted Wednesday that converts decimals to bits, and for mac users, a quick guide to getting started with Python will appear tomorrow.

Each of these bits is either 0 or 1, and goes towards storing a "double". From wikipedia.

We represent each (double precision) number x as (to use the lingo)

$x \approx sign*mantissa*2^{exponent-1023}$,

where the mantissa and exponent are each an integer.  We will deal with doubles, which means we get 64 bits, each either 0 or 1, to work with (see the nice figure above).  Then

1. The first bit (#63 above) records the sign of the number: 0 for positive, 1 for negative.
2. Bits 63-52 (11 bits) record an exponent E, the “size” of the number.
3. Bits 52-0 (52 bits) record the “mantissa” M of the number itself.  We always choose the mantissa to be as large as possible, which means the leading digit will always be 1, so we will actually squeeze out 53 bits of information by leaving out the leading 1.

Also notice that with 64 bits, we can hope to distinguish at most $2^{64}$ unique numbers.

Some nice hills in Montana.

As an example, we might store the number 1 as follows:
1. The sign is positive, so the first digit is 0.
2. E should be chosen so that the mantissa lies between $2^{53}$ and $2^{54}-1$.  Some inspection shows that E = 970 is the only choice, so the exponent is $2^{970-1023}= 2^{53}$. Note that in binary, 970 = 512+256+128+64+8+1 = 01111001001.
3. M should then be  $2^{53}$, stored in binary, with the leading 1 omitted.  This is just a string of 52 zeros.
4. So the 64 bits to store 1 is 0 01111001001 0000000000000000000000000000000000000000000000000000.

Now the problem: we cannot store most decimals exactly.  Let’s write the number 1.1.  Notice that if you type 1.1 into a Python prompt and hit return, 1.1000000000000001 is displayed, which is already a problem.  But anyways:

1. The sign is positive, so the first digit is 1.
2. Again we choose so the mantissa is in the correct range.  Again will be 970.
3. Now we want to choose M so that $M\cdot 2^{-51}$ is as close to 1.1 as possible, so we want to round the number $\frac{11}{10}*2^{51}$ to the nearest integer, and then encode that integer as a binary number.  This is easy to do with the Python command
 v,r = divmod(11*(2**51),10),
noting that the remainder is greater than 5 (namely, 6), so we should round up v = 9907919180215090 to v+1. Then
bin(v+1) outputs as a binary number.  In particular, we find that
1.1 “=” 0 01111001001 0001100110011001100110011001100110011001100110011001

This bench was mentioned in a sign in an earlier photo...

This is awful close, but not exactly it.  As a last two points:
1. The approximation is off by at most $2^{-53}$ as compared to the leading digit.  So the approximation to 1.1 was this close, whereas 4.1 might be wrong by $2^{-51}$, since the leading digit is about 4 times as big.
2. We can represent all integers between 0 and $2^{53}-1$ precisely.  After that, there is some error.

# An identity and floating point arithmetic

A geometric interpretation of arctan.

The other day I was looking at a string of calculations which were supposed to come out to zero.  Instead, the conclusion of these calculations was written as:

$\arctan{5} + \arctan{1/5}-\arctan{3}-\arctan{1/3}.$

Realizing that not only must this be wrong, but severely wrong, I plug it into a calculator just to make sure.  The results are:
1. Via google.com: 0.
2. Via wolframalpha.com: -5.551115123125783×10-17 is the answer given while typing it in, though after hitting “enter”, we get 0, but we also get a whole bunch of series expansions and continued fractions, which is odd.
3. Via MATLAB: -5.551115123125783e-17
4. Via Python: -5.5511151231257827e-17, which is just one extra digit on wolfram and MATLAB’s output.

First, the right answer is in fact 0.  Remember that arctan(x) takes the ratio of the opposite and adjacent sides of a right triangle (remember, side-angle-side completely determines a triangle), and returns the angle of the triangle, in radians.  Then it follows that since a triangle has pi radians total, and the right angle uses up pi/2 of those radians, we have the following identity for all nonzero x:

$\arctan{x} + \arctan{1/x} = \pi/2.$

It follows, then, that

$\arctan{5} + \arctan{1/5}-\arctan{3}-\arctan{1/3} = 0,$

though this is one of the most… interesting ways I have ever seen the number 0 written.

Tomorrow: more on why the calculators were wrong, and a rough discussion of exactly how computers calculate with floating point numbers.

# A little Python

Two filler photos for today. Neither of them contain pythons.

I had thought of starting today’s post with a picture of a snake, but I think yesterday featured enough unpleasantness.  We had two colloquia today, including a great talk by Luis Caffarelli on diffusion equations which had some great insights into regularity for these solutions.  But I’d like to talk about the second one, which was on the Frobenius coin problem, asking what the largest number  a certain number of coins could not make change for.  The idea is that if you lived in a (rather silly) country with only 3 and 5 cent pieces, you would be at a loss to give someone 7 cents.

But this is not the problem I want to talk about.  I want to think about countries that are moderately less silly.  They will always have coins to make any amount of change (this is equivalent to always having 1 cent pieces).  But the observation I seek to generalize is the following:

With [common] US coins, it is possible to have $1.19 in change, but not be able to make change for a dollar (i.e., 3 quarters, 4 dimes, and 4 pennies). Houston, looking positively bucolic. I want to imagine a country where a dollar is worth d pennies, and to investigate the worst choice of denominations for coins, so that a person can have a lot of money in coins, but not make change for a dollar. I also impose the following restrictions: 1. The country has at most 4 denominations of coin. 2. Each denomination of coin divides the value of a dollar. I wrote a short Python script to solve this. It turns out that as long as you have 4cent and 25cent pieces, you can hold$1.71 without being able to make change for a dollar.

The first time you get over $5 in coins is when 1 dollar = 210 cents, then if you have pennies, 30, 35, and 42 cent pieces, you can have$5.23 in coins without being able to break a dollar.  Similarly, if 1 dollar = 396 cents, you can have \$10.09 in coins, if coins have value 1,36,44 and 99.

I have put a plot below.  It reminds me vaguely of the plot of Euler’s totient function, and certainly a thorough description of a solution to this should have to do with how many factors the denominations of coins have, but I certainly do not have a rigorous explanation for why these two plots look similar.  The points along the bottom correspond to prime numbers, where condition 2 requires that we only have pennies.

A plot for value of a dollar, and the most change you can have without being able to break a dollar (assuming a really bad choice of denominations of coins).

Euler's phi function, for comparison.

# Writing down numbers

It has become a “go to” line for me, when discussing “what is math research?” with my (utterly fantastic) relatives or (pretty fantastic) non-mathematical friends that “someone’s gotta keep inventing new numbers.”  In an effort to mix a whole lot of influences for me, today’s post is inspired by three things:  The first the following dinosaur comic:

From Ryan North's fantastic Dinosaur Comics.

The second influence is from Edward Tufte, who loved (and I assume loves) using data to give micro/macro readings, like bar charts where you use numeric data instead of bars- the “bars” give a macro reading of how big the number is, while the number itself gives detail.

The 10th row of Pascal's triangle.

Finally, I really liked the wikipedia article on the binomial coefficient, in particular the graphic on the rows of Pascal’s triangle.  I’ve included in this post a printout of the 10th, 100th and 1000th row of the triangle.  For interest’s sake, the Python code I used to print these is

 def binom(n,k): if k > n: return 0 k = min(k, n-k) c = 1 for j in range(1,k+1): c = c*(n-(k-j))/j return c
 TOP = 1000 for k in range(TOP+1): print binom(TOP,k) 

100th row of Pascal's triangle. The numbers are written vertically, because captions look silly underneath thin figures.

I love that the simple act of writing down these numbers suggests that there is some pattern in the number of digits of the binomial coefficient.

The 1000th row. This actually took forever to make, since it is hard to take a picture of a text document. Turns out after all the compressing I did to get this, the "digits" are all just one or two pixels. Sigh...

For the curious, wikipedia reports that this is called “log concavity”, though a more satisfying answer (which I assume has already been found) might be to describe precisely and analytically the shape that the nth row of the triangle will make as n gets very large.  Likely you would first rescale everything to go from, say, -1 to 1, and the height of the kth row would be the log of nCk, rescaled by the same factor.  Anyways, how cool that 11 lines of code (or a dinosaur writing down 1000 integers) can raise such interesting questions!

# C++ and Python

Just some early thoughts as I got through the first 10 problems on Project Euler for the third time, using C++ this time:

• It is fast.  I knew that having a bit of code that would return a list of prime numbers was useful, and my first naive version of it can generate a list of 100,000 in just over 5 seconds, as opposed to ~11 with a more polished snippet in Python.
• I hate that everything happens at once.  When I fix programs in Python, I have been printing out information at various places to see if the code is doing what I expect.  Does not work in C++.
• The first few Project Euler questions have nice pencil and paper solutions, and it is sometimes hard to ignore those in order to practice programming.
So far, as long as everything compiles, I have not had to worry about time constraints at all, and don’t think I should.  I’ll optimize code later- now I’m just working out syntax.

Now I’m off to the Joint Meetings!