# Calculating the median in the St. Petersburg paradox

I am studying a recreational probability problem (which from the comments here I discovered it has a name and long history). One way to address the paradox created by the problem is to study the median value instead of the expected value. I want to calculate the median value exactly (not only find bounds or asymptotic values). I have found a certain approach and I am stuck in a specific step. I present my analysis and I would like some help on that specific step.

[Note: Other solutions to the general problem are welcome (however after the revelation of the long history I found a lot of material) but what I really want is to know the answer to the sub-problem that my approach raises.]

## The problem

We have the following game: I toss a coin as many times is needed to get tails. Then I count the number of consecutive heads that preceded (call it h) and I give you $2^h$ dollars. How much are you willing to pay to play such a game? In other words, what is the maximum buy-in for that game you are willing to pay? Note also, that we can play this game any amount of finite times (each time with you paying the buy-in).

One straightforward way to answer this is to calculate the expected value of one game. This should be the upper limit for the buy-in. The expected value is the infinite sum of the return of each case times the probability of each case. More specifically $$\sum_{i=0}^\infty (2^{i-1}\cdot\frac{1}{2^i}) = \sum_{i=0}^\infty \frac{1}{2} = \infty$$ This might seem counter-intuitive but it is true: Whatever constant and finite amount you bet per game, you are expected to win on the long run! Why is this so counter-intuitive though? Would you be willing to play this in practice with say 1000 dollars per game? The answer is no, because you would need an immensely large amount of games to actually win. So if we care about a more practical measure, the expected value is of no help. What we need is the median (or any other percentile value). If we know the median return for N games, we can at least know that if the buy-in is $\frac{median}{N}$, half of the possible cases you will lose and for half you will win. We will not know how much we will win or lose (we do have an upper bound on the losses though) but at least we know the chances to win or lose for a finite N number of games.

## Finding the median

So how do you calculate the median return from N games (or more generally any ith percentile)?

If we play only one game (N=1) then it is trivial. The median is 1. For N=2 it starts getting more complicated. With probability 0.25 we'll get back 1+1, with 0.125 1+2, with 0.125 2+1. These 3 cases already bring us to a total of 0.5, so the median is 3 (and so the maximum bet is 1.5 per game). For any N, how do we enumerate all the cases and find the 50% point (or any i% point)? I realized that this is (partly) an ordering problem. We do not want just to enumerate random cases, we have to order them, starting from the case with the smallest possible return, then getting the one(s) with the next smallest return and so on. As we are doing this ordering we are adding the probabilities of these cases. When we reach 50% (or i%) we stop. The return value for that case is our median value (ith percentile value). The ordering is where I am stuck.

### Sub-problem formulation

We can depict the possible space of returns with a matrix where the N columns are the N games and the infinite rows are the return for each game: $$\begin{array}{c} \text{row 1} \\ \text{row 2} \\ \text{row 3} \\ \vdots \\ \text{row i} \\ \vdots \end{array} \;\;\;\; \overbrace{\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ 2 & 2 & \cdots & 2 \\ 4 & 4 & \cdots & 4 \\ \vdots & \vdots & \ddots & \vdots \\ 2^{i-1} & 2^{i-1} & \cdots & 2^{i-1} \\ \vdots & \vdots & & \vdots \end{array}}^N$$

A series of N games consists of picking values for each column (i.e., picking a game outcome for each game). The smallest possible total return is when all game outcomes are 1. So total return = N. The next possible one is when we get one outcome from the second row (total return N+1). The next smallest total return is N+2 (2 game outcomes from the second row). Notice though that for total return N+3 we have two "configurations": 1) cases where we have N-3 outcomes from the first row and 3 from the second row, OR 2) cases where we have N-1 outcomes from the 1st row and 1 outcome from the 3rd row! So ordering is not such an easy process.

### Configurations vs. cases

Notice how I talked about "configurations" instead of individual cases. An individual case is a sequence of game outcomes (which are completely described by the game returns). For example a case of 4 games could be (1, 1, 16, 8) for a total return of 26. A configuration on the other hand is a more general construct which specifies how many outcomes we have from each row. A configuration completely determines the total return, but not the individual order that the outcomes happened. For example, the case given above is part of the configuration "2 outcomes from row 1, 1 outcome from row 4, 1 outcome from row 5". Cases (1,16,1,8) and (8,1,1,16) belong to the same configuration. From a configuration I can calculate how many distinct cases it has and what is the probability of each case. For example, for the configuration " $N_i$ outcomes from row i, $N_j$ from row j, $N_k$ from row k" we have:

The number of distinct cases is ${N\choose {N_i}}\cdot{{N-N_i}\choose{N_j}}\cdot{{N-N_i-N_j}\choose{N_k}}$

The probability for each of these cases is $2^{-(i\cdot N_i + j\cdot N_j + k\cdot N_k)}$

The total return value for any of these cases is $N_i \cdot 2^{i-1}+N_j \cdot 2^{j-1}+N_k \cdot 2^{k-1}$

The example above shows a configuration with 3 rows, just to get a taste of the complexity of the problem. I can generalise the formulas to find distinct cases, their probabilities and their total returns for any given configuration. The problem is ordering the configurations. Can we find an algorithm that orders and lists the configurations based on their total return value? Let's describe each configuration as a series of pairs {(x,i), (y,j), ...} where the first number of a pair denotes the row number and the second number of a pair denotes how many outcomes do we have from that row. For example, {(1,4), (3,1), (4,2)} means that we get 4 outcomes from row 1, 1 outcome from row 3, and 2 outcomes from row 4. This also means that we played 4 + 1 + 2 = 7 games. I manually computed the first terms of the ordered configurations list, for N games. I give the configuration(s) on the left and the total return on the right. Note that some total returns have more than one configurations that produce them.

$\begin{array}{ll} \text{Configurations} & \text{Total return} \\ \{(1,N)\} & N \\ \{(1,N-1),\; (2,1)\} & N+1 \\ \{(1,N-2),\; (2,2)\} & N+2 \\ \{(1,N-3),\; (2,3)\},\;\; \{(1,N-1),\; (3,1)\} & N+3 \\ \{(1,N-4),\; (2,4)\},\;\; \{(1,N-2),\; (2,1),\; (3,1)\} & N+4 \\ \{(1,N-5),\; (2,5)\},\;\; \{(1,N-3),\; (2,2),\; (3,1)\} & N+5 \\ \{(1,N-6),\; (2,6)\},\;\; \{(1,N-4),\; (2,3),\; (3,1)\},\;\; \{(1,N-2),\; (3,2)\} & N+6 \\ \end{array}$

If I can produce this order algorithmically then I will be able to calculate the median (or ith percentile) for any N.

I would also appreciate any help in formulating the problem in more accepted/mainstream terms. I believe that the formulation is valid and clear(?), but if we use a formulation from an established subfield maybe it will point to the solution too. Thanks!

-
Some more ideas can be found here: en.wikipedia.org/wiki/St._Petersburg_paradox. – Yuval Filmus Mar 20 '11 at 16:49
In fact if you play one game, then I think the median is between 1 and 2, since 1 has probability 0.5 of happening. – Henry Mar 20 '11 at 19:36

Feller was able to determine the median asymptotically, namely $N \log_2 N$, or $\log_2 N$ per game. He actually proved an even stronger concentration bound. This result can also be found in his "Introduction to Probability".

Economists are still coming up with new interpretations of this problem, known as the St. Petersburg Paradox, after the city in which Euler, who invented it, was working (the name might be due to Feller). See, for example, this article, which also suggests the same value.

-
Feller references earlier books, so the name is definitely not due to him. – Yuval Filmus Mar 20 '11 at 17:27
The Wikipedia article says "The paradox is named from Daniel Bernoulli's presentation of the problem and his solution, published in 1738 in the Commentaries of the Imperial Academy of Science of Saint Petersburg (Bernoulli 1738). However, the problem was invented by Daniel's cousin Nicolas Bernoulli who first stated it in a letter to Pierre Raymond de Montmort of 9 September 1713 (de Montmort 1713)." – Henry Mar 20 '11 at 19:17
You can find Feller's 1945 article at ProjectEuclid – Henry Mar 21 '11 at 7:17
Reading the linked articles, it seems that Feller had payouts twice those in this question (i.e. the shortest round would payout 2 rather than 1, the next shortest 4 rather than 2 etc.) so for this question the asymptotic medians would be half Feller's, which happens to be the $N \log_4 N$ or $\log_4 N$ suggestions which appeared in my simulations. That consistency is a relief as I had started to worry. – Henry Mar 21 '11 at 11:21
What is the stronger concentration bound you are mentioning? We know that $\frac {S_n}n - \log n\overset{d}\to G$ for an infinitely divisible $G$ with Pareto tail. Is this what you had in mind? – A.S. Dec 19 '15 at 4:46

It is easy enough to simulate this to get a reasonable estimate of the median values per game. For example using R

StPPmedian <- function(sampsize, numgames) {
median( rowMeans( matrix(2^floor(-log2(runif(numgames*sampsize))),
nrow=sampsize, ncol=numgames) ) )              }


you can get something like the following estimates of the median for 100,000 simulations each for various numbers of games:

> StPPmedian(100000,5)
[1] 2.4
> StPPmedian(100000,10)
[1] 2.9
> StPPmedian(100000,20)
[1] 3.4
> StPPmedian(100000,50)
[1] 4.08
> StPPmedian(100000,100)
[1] 4.59
> StPPmedian(100000,200)
[1] 5.115
> StPPmedian(100000,500)
[1] 5.78


and drawing these estimates of the median against the logarithm of the number of games certainly suggests some kind of logarithmic relationship is plausible, possibly with the median value per game close to $\log_4(N) + O(1)$ as also suggested in the following simulations

> StPPmedian(100000,4)
[1] 2.25
> StPPmedian(100000,16)
[1] 3.3125
> StPPmedian(100000,64)
[1] 4.296875
> StPPmedian(100000,256)
[1] 5.289062


In the comments I said that empirically the median appeared to be $\log_4(N)+O(1)$ and that to one decimal place the $O(1)$ term seems to be about $1.3$ for large $N$. A.S. suggested that the value of $O(1)$ depends on along which sub-sequence you take the limit. The following chart may demonstrate this: there seem to be visual patterns, combined with the noise resulting from simulation.

maxn <- 512 ; n <- 1:maxn ; meds <- rep(NA,maxn)
for (i in n){ meds[i] <- StPPmedian(1000000,i) }
plot(meds-log(n,4), ylim=c(1.2,1.4))


-
Thanks. I had done some simulation some time ago to get a feel, and the all other results point to a logarithmic relationship. But I am really interested to calculate the median exactly, using the ordering approach I described. – Thanassis Mar 21 '11 at 4:37
@Thanassis: It is probably very difficult to calculate the median exactly. At some point in mathematics, you realize that you have to make compromises sometimes. In the present case this means that you'll have to settle with some asymptotic expression. The first order statistics are given by Feller. You might ask what is the entire asymptotic series. This is the (mathematically) interesting question. – Yuval Filmus Mar 21 '11 at 6:28
While some mathematicians are interested in exactly calculating things, even in their case the results are not always explicit. In your case, for small $N$ one can probably calculate the exact value, but for moderately large $N$ this is probably hopeless. – Yuval Filmus Mar 21 '11 at 6:29
@Yuval I am not looking for a closed formula. I am looking for an algorithm that produces the ordering, which will in turn allow us to calculate the median (or any ith percentile). I do not think this is particularly difficult and maybe the problem is already solved in other settings. With this in place we can calculate medians for $N \approx 10^6$, something we cannot do efficiently with simulation. The asymptotic expressions are also interesting and I will study the material you suggested. – Thanassis Mar 21 '11 at 6:52
@Thanassis: I seriously doubt there is an algorithm which allows you to caclulate medians precisely for $N = 10^6$. We know the median sum will be over $10^7$ (I would guess the median average was around 11.2 or 11.3), so we have to consider prizes at least up to $2^{23}$. So I would guess that the number of possible patterns, adding up to the median sum or less, might more than the number of atoms in the universe, and each needs its own probability calculated. – Henry Mar 21 '11 at 12:44

You can estimate the median by considering events with probability more than $1/2$. For example, let $X_t$ be the event "exactly one game reached stage $t$". The probability that this happens is $$p_t = N 2^{-t} (1-2^{-t})^{N-1}.$$ For $t \approx \log_2 N$, we have $p_t \approx 1/2$. This means that with probability more than $1/2$, the total outcome is at least $2^t \approx N$. So the median is at least $\approx 2N$.

By considering the event "no game reached stage $t$", we get (in the same way) that the median is at most $N^2$. This estimate is probably very bad.

Both these values are not normalized - to get the value per round, divide by $N$.

Probably by devising better events you can get closer to the actual median. Your enumeration approach, on the other hand, seems to only work for small $N$, even if you did find a way to enumerate the probabilities in increasing order of total outcome. What I'm suggesting is very similar to this enumeration, only you don't have to be entirely accurate; the more accurate you are, the closer to the median you'll get.

-
This is a very interesting approach. Finding an approximate solution fast! I'll look more into it. But I think there is a mistake in your calculations/modelling. You say that a lower limit is 2N but in my example for N=2 I calculated the median to be 3. Did I do a mistake? I would still like to calculate the exact value. Ordering of the cases seems the only way. It should work for any N. I added a more advanced example which might illuminate the enumeration/ordering approach more. – Thanassis Mar 20 '11 at 15:48
My numbers are approximate - to get actual bounds you'll need to work out the calculations. Admittedly the bounds I get are not very interesting. – Yuval Filmus Mar 20 '11 at 16:46
@Thanassis: You are correct that with $N=2$ the probability of getting a total no more than 3 (or average per game of no more than 1.5) is 0.5, but it is therefore also true that the probability of getting a total no less than 4 (or average per game of no less than 2) is also 0.5, so the median is 3 or 4 or perhaps somewhere in between. – Henry Mar 21 '11 at 0:47
@Henry Yes, you are right. Because of the discrete nature of the total returns the median can be any number in an interval. For N=2 any number in [3, 4) can be the median. But we cannot include 4 because this is another distinct case that increments the cumulative probability by 0.125/2. I guess we can define median in such ways (using the equal or not) that we also accept the intervals (3,4) or [3,4]. – Thanassis Mar 21 '11 at 4:32

I had to write an ugly recursive algorithm in C to calculate the exact median for this game. Please find the first 30 exact median values in the graph at page 62 philpapers.org/archive/ERGTEO.pdf

-

The usual resolution to this paradox is to consider that whoever is paying you off cannot pay off an infinite amount of money. So you terminate your sum at the point you have all the bank's money. If that is \$1B, they can only pay 30 rounds, so the fair value is \$30. As it is logarithmic, it doesn't vary much as you change the estimate of what they can pay.

Added: You actually don't care about the order of results, just what they are. So in $N$ games, your chance of getting $N$ is just $2^{-N}$, your chance of getting $N+1$ is $N2^{-N}$, your chance of $N+2$ is $(N+N(N-1)/2)2^{-N}$ and so on. The generating function will keep track of this. If the cutoff is $2^m$, the generating function is $\left(\frac{x^{2^m}}{2^{m+1}} + \sum_{i=0}^m \frac{x^{2^i}}{2^{i+1}}\right)^N$ where the coefficient of $x^p$ gives the chance of winning $p$

-
OK let's assume that there is a finite amount to be won (1Billion, so up to 30 coin tosses, expected return is 30). This does not really answer the practical side of the question. Would you be willing to pay \$29 per game? You would have to play an very large amount of games to get your money back. What if you can only play 500 games? How much are you willing to pay then? How about for 1000 games? The median is an interesting measure and I would like to find it. – Thanassis Mar 20 '11 at 15:38
I do not need to know the order of the cases but I need to know the order of "configurations". Configurations was a structure that emerged after I realized that finding the cases for higher returns gets complicated. Initially I thought what you described but then I found the complication with return N+3 (explained in the question). I improved the structure of the question and how I explain the differences between cases and configurations. Have a look if you like and let me know if this helps. – Thanassis Mar 21 '11 at 4:46