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:
		count +=1
		coin = randint(0,1)
	return tot
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
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),:)]];
if N~=1
[A,B] = gasket(A,B,N-1);

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  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'
		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)
		d = 10**(len(n)-n.index('.')-1)
		n = int(d*N)
			d = 10**(len(n)-n.index('e')-1)
			n = int(d*N)
			d = 1
			n = N

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

	#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


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

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 Mine was in “colcarroll/mystuff/blog”. Once there, just type “python“, 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 0.
2. Via -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.