Toss Coin Until Head Appears Twice In a Row

Take a fair coin and toss it until a head appears twice in a row. There is no finite sample space for this experience, for it theoretically can go on infinitely. That said, the first few elements are: {H,H}, {T,H,H}, {T,T,H,H}, {H,T,H,H}. To determine the probability of tossing exactly 4 times, we can take the elements in the sample space that corresponds to only 4 tosses. In that case, they are {T,T,H,H}, {H,T,H,H}. The probability of tossing either combination is 1/16, and two such combination gives a total probability of 1/8.

To see the result simulated, the following codes on MATLAB can help to illustrate that with 1 million trials. Documentations are added to clarify the algorithm.

probability = 0.5;      %prob of tossing head
numTrials = 1000000;    %number of simulation trials
trialResults = zeros(1,numTrials);

for trial = 1:numTrials
    done = false;
    lastHead = false;
    toss = 0;
    while ~done         %this loop controls each trial
        roll = rand;
        toss = toss + 1;
        if roll < probability && lastHead
            trialResults(trial) = toss;
            done = true;    %each trial over when last trial was H, and this trial is also
        elseif roll < probability && ~lastHead
            lastHead = true;    %last trial was not H, but this trial is
        else
            lastHead = false;   %did not toss a H
        end
    end
end

sum(trialResults == 4) / numTrials

The results vary very little from 0.1250.

Showcase Showdown Analysis: Circular Reasoning and Oscillating Nash Equilibrium

The Showcase Showdown is a portion of the game-show The Price Is Right, in which the contestants spin the Big Wheel that has 20 sections randomly distributed, from 5 cents to $1.00 in 5-cent increments. The objective is to get as close to $1.00 as possible without going over, with one initial spin and an optional second spin. In the game, three contestants play the game to determine who has the highest value. What's the strategy in this game for the first contestant? An intuitive response may be to spin again for 50 or less on the first spin; keep the first spin if it's 55 or greater. Unfortunately, it's not as simple as that.

One strategy would be to spin the wheel the second time, only if the first spin resulted in less than the expected total outcome from the game, conditional to the outcome from the first spin. We can analyze the extreme discrete cases first to get a better understanding. If the first spin is 100, all second spins will make the total go over, resulting in total score of 0. In that case, the contestant will definitely keep the first spin of 100. On the other extreme, the expected total from the two spins, conditional to the first being 5, is 52.25. In that case, taking the second spin is better off. Here is the complete table, with the higher result reflecting the course of action pursued:

First Spin	2 Spins Avg Total	Higher Result
5	52.25	52.25
10	51.75	51.75
15	51	51
20	50	50
25	48.75	48.75
30	47.25	47.25
35	45.5	45.5
40	43.5	43.5
45	41.25	45
50	38.75	50
55	36	55
60	33	60
65	29.75	65
70	26.25	70
75	22.5	75
80	18.5	80
85	14.25	85
90	9.75	90
95	5	95
100	0	100

If the first spin were 40 or less, having a second spin will on average produce a better result. The average of the third column gives the expected result from the game: 63. However, there is one hole in this reasoning when applied to the game. This would work perfectly fine if the first contestant played the game for himself or herself, only concerned about maximizing the individual score given the risk-reward offset. Instead, of the three contestants, only the one with the highest result wins. If the first contestant got anything from 45 to 60, inclusive on the first spin, it would've been strategically better to keep it in the aforementioned reasoning. However, since the objective is to beat all other contestants, rather than maximizing individual scores, at that point of the game it may be more reasonable to spin again nevertheless.

This is where the circular reasoning kicks in. If the contestant's decision is to spin the wheel the second time if the first spin were less than 63, then the total expected value drops to 59.95. Essentially, while the contestant tries to base the individual decision given the overall expected result, the overall expected results depend exactly on the individual contestants' decision. It's a circular route of logic, and also illustrates game theory being applied. Having 40 and 45 be the cutoff would the best risk-reward optimization decision on the individual level. However when that is the baseline, the dominant strategy is then to use 63 as the decision's critical point. When players do that, they all incur more risk and distort the overall expected result downward. Would contestant then use 59.95 as the decision's critical point?

If the contestants did, the next critical point would be 61.3. Here, if 61.3 were the next critical point, it forces 59.95 to again be the overall expected value, just like 63 did. Therefore, the "Nash equilibrium" is a perpetual oscillation between 59.95 and 61.3. In the end, the only definitely conclusions are to spin again if the first spin is 40 or lower, and to keep the first spin if it's 65 or greater. Having 45 and 50 as the first spin stands in the grey area, and 60 is dead in no man's land, caught in between the oscillating equilibrium.

Maximum Possible "Jeopardy" Score and its Probability

What's the theoretical maximum score on the TV show Jeopardy? To get the maximum possible score, one has to ring first on all answers (or have a different first ringer answer incorrectly), answer everything correctly, bet the maximum amount on all occasions, and then a luck factor plays in. The Daily Double clues are located at the optimal spot. If all of the non-luck factors are obtained, what's the probability that yields all of the lucks in place to obtain the theoretical maximum score?

In the first round, there are 6 categories, with each with $200, $400, $600, $800, and $1000 clues. Add them up and multiply it by six: 6*(600*5) = 18,000. However, there is one Daily Double clue in this first round. Since the Daily Double clue annuls the monetary value originally associated with the clue, the Daily Double has to be hidden behind a $200 clue to minimize the annulling, and thus maximizing the overall score. Subtract 200 from 18,000 leaves 17,800. If everything is wagered on Daily Double and the answer is correct, that leaves 17,800*2 = $35,600 just from the first round.

In the Double Jeopardy round, everything is doubled: 6 categories, each with $400, $800, $1200, $1600, $2000. Furthermore, there are two Daily Doubles. Start with the $18,000 figure calculated from the previous round. Double that, since it's Double Jeopardy: 18,000*2 = 36,000. Now, we need to backtrack two of the $400 values to account for the Daily Double. That leaves 36,000 - 2*400 = $35,200 gained purely from Double Jeopardy, without the Daily Doubles. Don't forget to add the amount from first round, and that leaves 35,600 + 35,200 = $70,800. Now for the two Daily Doubles, as well as the Final Jeopardy, everything is wagered and answer correctly. That is 3 times of doubling the score: 70,800*2^3 = 70,800*8 = $566,400.

Now the question is, if suppose someone can ring first on all answers, answer everything correctly, and still bet the maximum amount on all occasions, what's the probability that the Daily Doubles will be located in the correct spot to allow this theoretical maximum score? Well, in the first round, it has to be hidden in one of the $200 clues. That's a 1/5 chance of that happening. Independently from that, the chance of the two Daily Doubles both hidden under $400 in Double Jeopardy is (6/30)*(5/29). The first one can be in any of the six $400 slots, out of the 30 overall clues. After that's taken care of, the second one has to be in any of the five remaining $400 slots, out of the 29 overall remaining clues. Therefore, the overall probability is (1/5)*(6/30)*(5/29) = 0.00689655172, which is exactly 1/145.

So overall, if someone can ring first on all answers, answer everything correctly, still bet the maximum amount on all occasions, and has the 1/145 luck chance, it is possible to get the theoretically maximum score of $566,400 on Jeopardy. To give a comparison, that value is over 1/5 of the amount Ken Jennings earned during his record-setting 74-win streak.

Sources:

• http://www.usatoday.com/life/people/2004-11-30-jennings-jeopardy_x.htm

4 Home-Run Game Probability

On Tuesday night, Texas Rangers outfielder Josh Hamilton hit 4 home-runs in one game. This was only the 16th time in Major League Baseball's history to see the accomplishment of such a feat, and the first since 2003. Hamilton also had a double in the game, going 5-for-5 overall and totaling 18 bases in a single game, which was only one shy of Major League record. The Rangers won the game 10-3 against Baltimore Orioles.

To best estimate the probability of this event, the following data for Hamilton during each of the past four seasons at the Rangers were retrieved:

Year	At Bat	Home Run
2008	624	32
2009	336	10
2010	518	32
2011	487	25
Sum	1965	99

Since Hamilton hit 99 home-runs in 1965 at-bats during the past four seasons, let's assume that the probability of hitting a home-run is 99/1965 = 5.038%. On Wednesday night, Hamilton had 5 at-bats, and hit 4 home-runs. The probability of hitting 4 home-runs in a game with 5 at-bats will be the product of the following terms:

• nCr(5,4) to indicate the number of combinations
• (5.038%) ^ 4 to indicate the 4 home-runs
• (1-5.038%) to indicate the 1 non-HR at-bat

Using this MATLAB line, the product is calculated to be nchoosek(5,4)*(99/1965)^4*(1-99/1965) = 3.0592e-005. That is around 1/32,688. Let's not forget the double he had. If we factor that into the calculation as well, the probability turns out to be 1/487,968.

To put in perspective, each season is 162 games. During each game across the league, there are 30 teams and 9 players who come to bat. That is 30*162*9 = 43,740 total player-performances in a season, assuming uniformity. Already that number is greater than 32,688. This means that if every player in MLB had the HR-hitting probability as Hamilton, the league should expect to see occurrence like Tuesday night's every season. However, clearly that isn't the case. Four-HR games are more rarer than perfect games or no-hitters. In baseball, only 20-K games and unassisted triple players have been less frequent.

Sources:

• http://espn.go.com/dallas/mlb/story/_/id/7907243/josh-hamilton-texas-rangers-hits-four-2-run-hrs-baltimore-orioles
• http://espn.go.com/mlb/player/stats/_/id/4652/josh-hamilton

Calculating Type II (β) Error

In hypothesis test, β is the probability of failing to reject the null hypothesis when it is actually false. Let's first consider the two-sided case to see how it is calculated. Begin with the following given conditions:

• Let x̄(n) = 174.5, s(n) = 6.9, n = 50, α = 0.05
• H0: μ = 175, H1: μ ≠ 175
• Compare against alternative hypothesis μ = 173

At α = 0.05, we are concerned with z = 1.96 and z = -1.96. Turn these z values into critical values for x̄:

• Recall that for samples (for both two- or one-sided cases), z = (x̄-μ) / (s(n) / sqrt(n))
• 1.96 = (x̄ - 175) / (6.9 / sqrt(50)) --> x̄ = 176.913
• -1.96 = (x̄ - 175) / (6.9 / sqrt(50)) --> x̄ = 173.087

Under standard conditions, the null hypothesis would be accepted if x̄ is between those values. But what if the mean is really 173? That is the alternative hypothesis we are considering. To do that, use a normal distribution with 173 as the mean. We then need to find the z-values for 173.087 and 176.913. To do that:

• (173.087 - 173) / (6.9 / sqrt(50)) = 0.0892
• (176.913 - 173) / (6.9 / sqrt(50)) = 4.0100

From here, β is just the area under the normal distribution curve between the two critical values. Using MATLAB, normcdf(4.0100)-normcdf(0.0892) = β = 0.4644. The power of the test, or probability of correctly rejecting the null hypothesis when it is false, is 1-β = 0.5356.

What about an one-sided case? It is similarly computed. Consider the following given conditions:

• Let p̂ = 0.535, n = 1000, α = 0.05
• H0: p = 0.50, H1: p > 0.50
• Compare against the alternative hypothesis p = 0.52

Remember that for the one-sided case, α = 0.05 gives the critical value of z = 1.645. Again turn the value into a value for p̂, similar to the way we solved for x̄ in the example above :

• Recall that for proportions (for both two- or one-sided cases), z = (p̂-p) / sqrt(p*(1-p)/n)
• 1.645 = (p̂-0.5) / (sqrt(0.5^2/1000)
• p̂ = 0.52601

Under standard conditions, the null hypothesis would be accepted if p̂ is 0.52601 or less. Now use a normal distribution with 0.52 as the mean. We want the area to the left of 0.52601, because this is the area that will not be rejected. Now, carefully note that if our H1 were p < 0.50, we would be concerned with the area to the right, instead of left, of 0.52601. Back in the scenario here, we now need to turn 0.52601 to a z-score.

• Do not forget that p = 0.52 now, instead of 0.50
• (0.52601 - 0.52) / sqrt(0.52*0.48 / 1000)
• = 0.3804

To obtain the area to the left of that z-score, simply run normcdf(0.3804) = 0.6482 = β. The power, or the probability of correctly rejecting the null hypothesis when it is false, would again be 1- β = 0.3518.

Probabilistic Payoff from Rushing to Train Stations

See related post: Expected Wait Time given Multiple Trains

Is it worth it to run and rush to the train station? It might payoff if running there gets us on an earlier train. But it might not payoff, as we simply end up waiting a longer time for the next train. Does it make a difference what time of the day it is?

To answer this question, we need just two parameters: the frequency of the arrival event, and the extent to which the rushing shortens the time getting onto the site. Let's use the train example for easier understanding. Suppose that the train comes every 10 minutes. As we saw on the previous post, the average wait time until the first train is 5 minutes. Now, suppose that we run to the train station and get there 3 minute earlier than if we had walked.

It turns out that this answer is much simpler than the calculations from the last post. The average wait time is still 5 minutes. However, the probability that we end up catching an earlier train now is 3/10. We can easily see that we'll be able to catch an earlier train, if a train comes within those first 3 minute that we wouldn't have been there otherwise. Given 10-minute frequencies, that probability is simply 30%. This helps us answer the question of whether or not rushing to the arrival site will benefit. Using the variables defined earlier, the probability of benefit is simply (shortened time) / (frequency). In the event that (shortened time) > (frequency), that probability simply caps at 100% certainty.

That simply results also helps us answer whether time of the day affects this decision. Since during rush-hours the trains come more frequent (lower frequencies), rushing to the station will be more likely to be beneficial then. It makes intuitive sense, because in the same 3 minutes or whatever the shortened time may be, there will be a higher chance that a train comes during that period then.

Expected Payoff from Mega Millions

Lottery games usually advertise the enormous amount of money that can be won, and the luxurious life that such winning would bring. But everyone knows that the chance of winning the jackpot is miniscule. So in this battle of extremities, does the large amount of money or the small chance win out?

Information taken right from the Mega Millions website reveals that it takes $1 to play a game. The calculations put aside for each of these possible prizes, here are the monetary amounts of the prizes with their associated chances of winning:

Prize	Chance ^ (-1)
Jackpot	175,711,536
$250,000	3,904,701
$10,000	689,065
$150	15,313
$150	13,781
$7	306
$10	844
$3	141
$2	75

The ultimate jackpot prize value varies consistently, but at present is estimated at $65 million. To get the expected payoff, simply multiply each earning by its likelihood, and sum the products. This can be easily done in Excel.

When the jackpot is $65 million, the expected payoff is only 55.18 cents. In fact, it takes slightly less than $144 million for the expected payoff to just reach $1. The jackpot recently reached record high of $656 million. With that, the expected payoff was slightly over $3.91. However, tax was not considered in any of these situation. Furthermore, as in the case of the recent record high, the winning amounts are split if more than one winner arises. In all, the small probability of winning seems to dwarf the large amounts of potential winning.

Sources:

• http://www.megamillions.com/howto/
• http://abcnews.go.com/US/wireStory/ill-mega-millions-winner-claim-jackpot-share-16160742#.T45qo9XZYS4

Expected Wait Time given Multiple Trains

Suppose a train comes every 10 minute, and we go to the train station at a random time. How much time do we expect to wait before the train comes? In this scenario, it's 5 minutes, because the probability is uniformly distributed between 0 and 10 minutes. This situation highly resembles waiting at the 116th Street - Columbia University station for a Downtown 1 train. Sometimes we get lucky and don't wait much; other times, we luck out.

Now suppose we take the train down to 96th Street and want to transfer to the 2 or 3 train. Suppose each 2 train and each 3 train comes every 10 minutes, independent of each other (this is a bit unrealistic, since they can't both come at the same moment on the same tracks). Anyhow, say if our destination is Chambers Street, we don't care whether we take the 2 or the 3 train; we just want to get on the train that comes first. What's the expected wait time here?

The expected time here is the time that elapses when neither the 2 train nor the 3 train comes. Let's backtrack first:

• Let f(x) be the probability density function of waiting x minutes at the platform. The function f(x) is simply 1/10 for x = [0,10], and 0 for the rest of the values of x.
• Let F(x) be the cumulative density function, which is simply the antiderivative of f(x) at F(x) = x/10 for x = [0,10]. 
• In example, F(7) = 70% means that probability that waiting time is 7 min or fewer is 70%. Or in another word, 1-F(7) = 30% means that there's a 30% chance that waiting time is 7 min or greater.

Let G(x) be the waiting time when there are two trains, each with f(x) = 1/10 for x = [0,10]. Then, the term 1-G(x) denotes the probability that neither train comes within x minutes. That happens with neither of the trains comes within x minutes. Since the two trains come independently in this exercise, 1-G(x) = (1-F(x))*(1-F(x)). Using the example of x=7 again, 1-G(x) = (30%)^2, so G(x) = 1-.3^2 = 91%. This means that there's a 91% chance that one of the trains will come within 7 minutes. So we saw that G(x) = 1-(1-F(x))^2. Recall that F(x) = x/10 for x = [0,10]. Therefore, we have G(x) = 1-(1-x/10)^2. Here's the graph of that:

What's the average waiting time? We need the probability density function g(x), which is just the derivative of G(x), and then integrate x*g(x) over the variable x from 0 to 10. Upon executing that command on MATLAB:

syms x;
y = 1-(1-(x/10))^2;
pdf = diff(y,x);
average = int(x*pdf,x,0,10)

The average time comes out to be 10/3 minutes. To verify that, a simulation of 100 million trials was run on MATLAB. The code isn't shown here (maybe shown in a future post), but the result came dead on at 3.3334.

What if there are more than 2 trains? With a frequency that isn't every 10 minutes? No problem. Look back at G(x) = 1-(1-x/10)^2. If we have understood what each number and variable has stood for, we should have no problem deducing that in the general situation with n trains, each independently coming every f minutes, the cumulative density function is G(x) = 1-(1-x/f)^n. To get the average, differentiate G(x) by x to get g(x), and integrate x*g(x) over x from 0 to n. Just for sake of curiosity, suppose it's rush hour at Chambers Street, and we can catch either one of the 1/2/3 trains to go back Uptown. If each train's frequency is 5 minutes, the average expected waiting time, following the formulas investigated here, is only 1.026 minutes.

For reference, the source of inspiration for the investigation of this subject was the countless amount of time spent waiting on the Subway platforms.

Binomial vs. Poisson vs. Normal Distribution

Consider this scenario: Given 100 bottles, each of which independently has 10% chance of being defected, what's the chance that up to 0 will be defected? Up to 1 will be defected? Up to 2? 3? ... Change to a general case, given (n) samples, each of which independently has probability (p) of returning true, what's the probability that up to (k) of the samples will return true? The binomial distribution gives an exact answer, while Poisson and normal distribution help to approximate an answer, with varying accuracy depending on the scenario.

Write MATLAB code as follows. In this case, n = 100, p = 0.05, and k = 0 to 10, but all of those values can be easily changed manually. The first column of the vector were manually-chosen (k) values for which the cumulative probability values want to be calculated:

result = 0;
table = zeros(11,4);
table(:,1) = [0;1;2;3;4;5;6;7;8;9;10];
n = 100;
p = 0.05;

for i=1:11
    table(i,2) = binocdf(i,n,p);
    table(i,3) = poisscdf(i,n*p);
    table(i,4) = normcdf((i+0.5-n*p)/(n*p*(1-p))^0.5);
end
table

Now use and tweak the program (value of n, p, and table(:,1)) above to run through two cases. In case 1, n is big while p is small. In case 2, p is relatively large. The combination of MATLAB results and Excel calculations were used to produce this table:

Case 1: n = 200, p = 0.02, λ = 4

k	Binomial	Poisson	Normal	Pois Error	Norm Error
0	0.0176	0.0183	0.0385	3.98%	118.75%
1	0.0894	0.0916	0.1034	2.46%	15.66%
2	0.2351	0.2381	0.2243	1.28%	-4.59%
3	0.4315	0.4335	0.4003	0.46%	-7.23%
4	0.6288	0.6288	0.5997	0.00%	-4.63%
5	0.7867	0.7851	0.7757	-0.20%	-1.40%
6	0.8914	0.8893	0.8966	-0.24%	0.58%
7	0.9507	0.9489	0.9615	-0.19%	1.14%
8	0.9798	0.9786	0.9885	-0.12%	0.89%
9	0.9925	0.9919	0.9973	-0.06%	0.48%
10	0.9975	0.9972	0.9995	-0.03%	0.20%

Case 2: n = 100, p = 0.4, λ = 40

k	Binomial	Poisson	Normal	Pois Error	Norm Error
15	0	0	0	#DIV/0!	#DIV/0!
20	0	0.0004	0	#DIV/0!	#DIV/0!
25	0.0012	0.0076	0.0015	533.33%	25.00%
30	0.0248	0.0617	0.0262	148.79%	5.65%
35	0.1795	0.2424	0.1792	35.04%	-0.17%
40	0.5433	0.5419	0.5406	-0.26%	-0.50%
45	0.8689	0.8097	0.8692	-6.81%	0.03%
50	0.9832	0.9474	0.984	-3.64%	0.08%
55	0.9991	0.9903	0.9992	-0.88%	0.01%
60	1	0.9988	1	-0.12%	0.00%
65	1	0.9999	1	-0.01%	0.00%

These data show that Poisson distribution is a better approximation when p is small, while normal distribution is a better approximation when p is large. As the numbers here were copied from MATLAB onto Excel, rounding errors have be distorted the percentage error calculations a bit.

Mixed Strategy Game

Two players A and B play a game, in which both has a red and a blue marble. They simultaneously present a marble to each other. If both present red, A wins $3. If both present blue, A wins $1. Otherwise, if the colors don't match, B wins $2. Who is in the better situation here?

		A
		Red	Blue
B	Red	A: $3	B: $2
	Blue	B: $2	A: $1

If both players choose red or blue with equal probability, the expected payoff for both players is identical at $1. It's true that Player A has a greater variance in the results, and therefore Player B is better for risk-averse individuals. However, let p be the probability that the other player chooses red. Here in the mixed strategies:

• 3p = (1-p) --> A would be indifferent about the decision if B had 1/4 chance of choosing red
• 2p = 2*(1-p) --> B would be indifferent about the decision if A had 1/2 chance of choosing red

What does this mean? Although there is no dominant strategy, Player A clearly prefers playing red (expected payoff = 1.5, against 0.5 for blue), if Player B chooses randomly. However, knowing this, Player B would be more likely to choose blue. 3/4 chance of choosing blue for Player B finally would make Player A indifferent. What about the expected payoff in this mixed strategy?

• Player A: (1/4)(1/2)(3) + (3/4)(1/2)(1) = 3/4
• Player B: (1/2)(1/2)(2) + (1/2)(1/2)(2) = 1

So in this mixed strategy, whereby Player B recognizes that Player A prefers playing red, Player B plays blue more frequently. Although Player A had preferred red under random decisions from B, now Player A is better off playing blue if B also puts blue.

Finally, to simulate the game on MATLAB, this script was written. This case simulates the mixed strategy where B has 25% chance of choosing red. In the original consideration of random selection, simply change the value of the variable RedprobB.

trial = 1;
numTrials = 1000000;
payoffA = zeros(1,numTrials);
payoffB = zeros(1,numTrials);
RedprobA = 0.5;
RedprobB = 0.25;
while trial <= numTrials
    testA = rand;
    testB = rand;
    if(testA < RedprobA & testB < RedprobB)
        payoffA(trial) = 3;
    elseif(testA >= RedprobA & testB >= RedprobB)
        payoffA(trial) = 1;
    else
        payoffB(trial) = 2;
    end
    trial = trial + 1;
end

sumA = sum(payoffA);
sumB = sum(payoffB);
disp(sumA)
disp(sumB)

If both players choosing red and blue with 50% chance, here are some of the results of the sum of the payoff over 1,000,000 trials:

• A: 999,425; B: 1,000,670
• A: 999,662; B: 1,000,664
• A: 1,001,261; B: 999,018

The results closely aligns with the expected payoff per trial of 1. Now, change back the code so that Player B has 25% chance of choosing red. Here are some of the results:

• A: 749,312; B: 1,001,864
• A: 751,635; B: 997,710
• A: 749,364; B: 999,236

Again, the results closely aligns with the revised expected payoff per trial of 0.75 for Player A and 1 for Player B. Now, taking a step further, if Player A knows that Player B will play red 1/4 of the time, what can Player A do? By changing the value of the variable RedprobA, here are some observations:

• RedprobA = 0.01: A earns about 0.75 per trial, B about 0.50
• RedprobA = 0.25: both A and B earn about 0.75
• RedprobA = 0.75: A earns about 0.75 per trial, B about 1.25
• RedprobA = 0.99: A earns about 0.75 per trial, B about 1.50

So it looks like if Player B is fixed at 1/4 probability of picking red, Player A can't do much to improve its own expected payoff. This shouldn't come at a surprise, since back earlier when p=1/2 for A and p=1/4 for B were calculated, that was the Nash equilibrium. Neither player has an incentive to switch strategies, given that the other won't. However, Player A here is able to affect the payoff of B, but it was assumed for this problem that the payoff of the other player doesn't factor into any decisions.

Probability of Consecutive Coin Flips

Given (n) coin flips, what's the probability of getting at least one pair of consecutive heads?

If n = 2, the probability is 1/4. If n = 3, the probability is 3/8 (HHH, HHT, THH). If n = 4, the probability turns out to be 8/16. The denominator of the probability fraction, in its unsimplified form, will be 2^n. The challenge is to find the numerator. For simplify in wording, here are the phrase abbreviations:

• "at least a pair of consecutive heads" --> "HH" 
• "no pair of consecutive heads" --> "!HH"
• "heads" --> "H" and "tails" --> "T"

Let p(n) be the probability that first (n) flips feature HH. This is what we're ultimately looking for. As with many cases in probability, it may be easier to define and use the complimentary probability q(n) as 1-p(n). Moreover, let's define q'(n) as simply the denominator of the unsimplified fraction q(n). So q'(2) = 3, q'(3) = 5, q'(4) = 8, and so on. To get started, take the first values of (n) and list out the possibilities of combinations featuring !HH.

n = 2	n = 3	n = 4
HT	HTT	HTTT
TH	THT	THTT
TT	TTT	TTTT
	TTH	TTHT
	HTH	HTHT
		TTTH
		THTH
		HTTH

While the ordering wasn't particularly important at this stage, notice that the first three rows of n=3 was conveniently laid out so that they are simply addition of one more T from the same row under n=2; same for going from n=3 to n=4. Now take a looks specifically at n=3. The first toss can be H or T. For 3 flips featuring !HH:

• If the 1st flip is T, there are 3 possible combinations (THT, TTT, TTH)
• If the 1st flip is H, there are 2 possible combinations (HTT, HTH)

The combinations of those 5 combinations happen to be the 5 combinations listed in the column above for n=3. Now start from n=2:

• For 2 of those combinations ending with T, each will produce 3 further possible combinations in 3rd and 4th flip to keep the !HH property
• For the remaining combination ending with H, it will produce 2 further possible combinations instead

Add 2*3 + 1*2 = 8, and it is no surprise that q'(4) = 8 as established. How the question remains: how did we just know that for n=2, two combinations ended in T?

The answer actually indirectly lies in q'(3). Observe the column of combinations above for n=2 and n=3. In particular, look at the last two rows. Similarly, compare n=3 and n=4 rows and observe the last three rows of n=4. What are those "extra" rows without something to the left of it in the table? They all end in H, and the rest of the term (without the final H) are exactly the combinations that end in T for the previous column.

Let's formulate some ideas for q'(5). Take q'(3) = 5 to start.

• Each of those 5 combinations will produce at least two possible combinations (of 4th and 5th flip) to keep the !HH property ongoing at n=5
• Furthermore, some portion of those 5 combinations will produce a third combination, if it ends in T

How did we how many ends in T? It's q'(4) - q'(3), or the number of those "extra" rows.

• Summing them together, q'(5) = 2*q'(3) + (q'(4) - q'(3))
• That simplifies to q'(5) = q'(3) + q'(4)

Suddenly, it looks very familiar. That is exactly the Fibonacci sequence that q'(n) is following, with q'(1) = 2, q'(2) = 3, etc. Precisely, given the initial offset, the Fibonnaci sequence F(n) = q'(n-2) or F(n+2) = q'(n). But recall, we're trying to find p(n), and in the process, defined q(n) and q'(n). To convert q'(n) = q'(n-1) + q'(n-2), recall that q'(n) is simply numerator of q(n), whereby the denominator is 2^n.

• Multiple each term by either 2^n, 2^(n-1), or 2^(n-2)
• That ends up being 4*q(n) = 2*q(n-1) + q(n-2)
• Simplify to q(n) = (1/2)*q(n-1) + (1/4)*q(n-2) with conditions that q(1) = 1, q(2) = 3/4
• Finally, p(n) = 1 - (1/2)*(1-p(n-1)) - (1/4)*(1-p(n-2))

Alternatively, it may just be more convenient to leave the result in Fibonacci sequence, since those numbers are more easily accessible. In that case, it goes from:

• Given F(n+2) = q'(n)
• q(n) = F(n+2) / (2^n)
• p(n) = 1 - F(n+2) / (2^n)

In either case, q(n) gives the probability of featuring no consecutive heads in the first (n) flips; p(n) gives the probability of having consecutive heads. To give some numerical sense, Java coding was written for the modified Fibonacci method. Beware that given the recursive method, the program will run quite slowly as (n) gets around 40.

  int until = 0;
  System.out.print("n = ? ");
  until = reader.nextInt();
  for(int i=1; i<=until; i++)
  {
      System.out.println(i + ": " + (1-fib(i)));
  }

  public static double fib(int run)
  {
      if(run == 1)
          return 1;
      else if(run == 2)
          return 0.75;
      else
          return (0.5)*fib(run-1) + (0.25)*fib(run-2);
  }

The condensed results are as follows:

n	p(n)
1	0.00%
2	25.00%
3	37.50%
4	50.00%
5	59.38%
6	67.19%
7	73.44%
8	78.52%
9	82.62%
10	85.94%
15	95.13%
20	98.31%
25	99.41%
30	99.80%
35	99.93%
40	99.98%
45	99.99%

As expected, the probability quickly climbs up as the number of flips increase, and asymptotically approaches certainty.

Variation of St. Petersburgh Paradox Simulation

Related previous posts: Simulation of St. Petersburg Paradox

In the traditional St. Petersburg Paradox, there is 1/(2^n) chance of an outcome of (2^n), for all positive integers n, summing up to an expected value of infinity, but with realistic outcomes far lower than that. What if the base is changed to an integer other than two? In this exercise, bases of 2, 3, 4, 10, and 100 were used to run similar simulations, each with 100,000 trials. To perform this variation, only two lines of codes needed to be changed from the original Java code. For base 3, for example:

• "while(flip < 0.5)" becomes "while(flip < 1.0/3)"
• "double result =  Math.pow(2,times)" becomes "double result =  Math.pow(3,times)"

After running the simulation, the results and their analysis as follows:

Base	Average/Base	Max/Base	St. Dev	Log(max)/Log(base)
2	9.75	65536	784.25	17
3	7.49	59049	724.34	11
4	6.47	65536	956.40	9
10	5.48	100000	3315.75	6
100	3.46	10000	12282.46	3

Given the different bases, the average result / base value becomes the most significant output to compare. As the value of the base increases, the normalized average decreases, since the probability of getting the lowest value dramatically increases. For base of 100, the result will simply be 100 for 99% of the time. While it's true that the output will also be dramatically greater, the clustering around high-probability, low-value region outweighs the rarity of low-probability, high-value regions.

The max result / base value is a bit harder to decipher. This is the output value were the game started with output of 1, rather than the value of the base, offering a normalized maximum output value. There doesn't seem to be a clear-cut trend of these values as the bases change. The lack of trend is reasonable, given that it only takes one rarity event to record the maximum. The standard deviation drastically increases as the base increases, a more obvious reflection of the greater variance in the output as the base increases. The logarithm of maximum value / logarithm of base is equivalent to the maximum number of times that the game ran for. The decreasing trend is reasonable, given that the probability of extending the game decreases dramatically as the base increases (50% for base 2, 1% for base 100).

The last three trends (or two trends plus the lack of trend for max/base) are reasonable. But the average result / base value is still the most intriguing, since the expected values, regardless of the base, are still all infinity. Instead, as expected, the greater the variance (bigger bases), the further the deviation from the expected value to the realistic output value (smaller normalized results).

Simulation of St. Petersburg Paradox

St. Petersburg Paradox exemplifies a situation where the expected value of an outcome doesn't reflect the realistic outcomes. The game can be played by flipping a coin until heads comes up. Let the outcome be 2 to the power of the number of coin flips. The derivation won't be shown here in detail (but is available on the link below), but essentially, there is 1/(2^n) chance of an outcome of (2^n), for all positive integers n. The expected payoff for each integer n is therefore 1, and the sum of the expected payoffs, the expected value of this game, is therefore infinity.

Realistic outcomes are far from this expected value. One can briefly visualize this by realizing that there's 50% chance of an outcome of 2, 25% chance of an outcome of 4, 12.5% chance of an outcome of 8, etc. Simply put, the chances of very high outcomes are very small. Let's use programming to simulate numerous rounds of this game. The Java code for the simulation method is as follows:

public static int simulation()
 {
     int times = 0;
     double flip = 0;
     while(flip < 0.5) //reflecting the 50% chance of getting heads
     {
         flip = Math.random();
         times++;
     }
     double result =  Math.pow(2,times);
     return (int)result;
 }

Writing a for loop code, this method was called 1,000,000 times, with the output values exported onto Microsoft Excel. The values were sorted and counted, and the results are as follows, along with the expected frequency of the different outcomes.

n	Outcome (2^n)	Frequency	Expected Freq	Deviation
1	2	500,592	500,000.00	0.118%
2	4	249,969	250,000.00	-0.012%
3	8	124,870	125,000.00	-0.104%
4	16	62,262	62,500.00	-0.381%
5	32	30,951	31,250.00	-0.957%
6	64	15,741	15,625.00	0.742%
7	128	7,843	7,812.50	0.390%
8	256	3,883	3,906.25	-0.595%
9	512	1,900	1,953.13	-2.720%
10	1,024	993	976.56	1.683%
11	2,048	484	488.28	-0.877%
12	4,096	233	244.14	-4.563%
13	8,192	144	122.07	17.965%
14	16,384	71	61.04	16.326%
15	32,768	34	30.52	11.411%
16	65,536	12	15.26	-21.357%
17	131,072	11	7.63	44.179%
18	262,144	3	3.81	-21.357%
19	524,288	4	1.91	109.715%

Even with 1 million trials, the greatest output value was only 524,288, which reflects the rare instance of 18 straight flips of tail before finally getting a heads. This only happened 4 times out of the 1 million trials, and was even considered strong positive deviation from the expected frequency of 1.91. The expected value of this game (expected average value) is still infinity, as the math demonstrates.

In this simulation of 1,000,000 trials, the average output value was 20.496084. The results are, as expected, very skewed to the right. If we let the output value measure monetary amount, it's interesting to note that the top 996 trials (0.0996%) contains 51.30% of the values. The top 0.7772% contains 65.86%. The bottom 96.8644% contains only 24.33%.

Sources:

• http://plato.stanford.edu/entries/paradox-stpetersburg/

Probability Limit

At a Science Café talk about mammals tonight, I heard this statement: "even if the chance of it happening is one in one-million in a given year, the likelihood is high given a span of one-million years." Is that really true? What about in the general case of "if the chance of an event happening is 1 in n, what's the likelihood that the event will occur at least once for n number of trials?"

The expression for that probability is 1 minus the chance that the event will not occur at all during the n number of trials. The chance of the event not occurring each time is (n-1)/n. Raise that to the n-th power for the number of iterations, and subtract from 1, which represents certainty, and the probability that the event happens at least once is:

We wish to find the limit as n approaches infinity. Working out the L'Hopital's rule by hand can be somewhat cumbersome with the exponential term, but with the help of Mathematica, we determine the limit to be:

The answer approximates to = 0.632120559. As n increases, the output probability decreases. However, n only needs to be greater than 485 for the output probability to match the first 3 decimal places after rounding. Of course for any given probability, as the number of trials increase indefinitely, the overall probability of the event occurring at least once approaches certainty. But in the special circumstance here, where the number of trials (n) match the inverse of the probability of each event (1/n), the overall probability approaches (e-1)/e.