Great problem solved - original material

[BruceZ]

A couple of months ago, someone asked for the probability of getting 3 AA in a row in 250 hands. This gave rise to the more general question which was phrased as:

Do you have the formula for "at least 1 streak of at least m occurences in n trials"? I don't. Show us yours.

This is an important problem in the theory of runs. At the time, I had given a number of approximations which were very close to the exact answer. Last year, I discovered an exact solution to the particular problem for m=2 and probability of success = ½ in terms of the Fibonacci series. Since that time I have had a project on the back burner to extend this exact solution to the general question asked above. I have now completed this project, and I have not one but two different exact closed form solutions to this problem, and I have verified unequivocally that they are both correct. Although accurate approximate solutions can be found by elementary methods, the exact solution is of a very different nature than one may initially expect, and the solution is quite elegant and beautiful. Hence I am now prepared to “show you mine” as requested.

To restate the problem, we want the probability P(m,n,p) of at least 1 occurrence of a run of of length at least m of an event with probability p in n trials. For the question of 3 AA in a row in 250 hands, m = 3, n=250, p=1/221. For a coin flip problem, p=1/2.

Exact closed form solution 1

The probability P(m,n,p) of at least 1 occurrence of a run of length at least m of an event with probability p in n trials is given exactly by:

P(m,n,p) = Sum{k = 1 to n} F[k]

Where F[k] is given by the recurrence relation:

F[k] = F[k-1] - (1/p - 1)*F[k-m-1]*p^(m+1) for k > m+1

And the initial conditions:

F[k] = 0 for k< m
F[m] = p^m
F[m+1] = (1-p)*p^m

The F[k] can easily be computed and summed in Excel. For the 3 AA problem simply set:

A1 = 0
A2 = 0
A3 = 1/221^3
A4 = (220/221)*(1/221)^3
A5 = A4-220*A1/221^4
.
.
.
Copy the last formula down to cell A250. Then compute =sum(a1:a250) for the exact probability. An Excel spreadsheet is available upon request.

For m=2 and p=1/2, this reduces to

P(2,n,1/2) = sum{k = 2 to n} F[k-1]/2^k

Where F[k] is the Fibonacci series as derived previously, and this simplifies too:

P(2,n,1/2) = 1 - F[n+2]/2^n

Subsequent to this derivation, I have found a relatively recent reference to the Fibonacci series solution for this special case in the literature; however, I have not found the closed form exact solution to the more general problem which I have given above . There is a slight chance that I am the only one in possession of it.


Exact closed form solution 2

The probability P(m,n,p) of at least 1 occurrence of a run of length at least m of an event with probability p in n trials is given exactly by:

P(m,n,p) = v[m+1]

Where v is an m+1 dimensional vector given by the m+1 x m+1 matrix product:

<pre><font class="small">code:</font><hr>
[ 1-p 1-p 1-p . 1-p 0 ] ^n [ 1 ]
v = | p 0 0 . 0 0 | | 0 |
| 0 p 0 . 0 0 | * | 0 |
| 0 0 p . 0 0 | | 0 |
| … . … | | … |
[ 0 0 0 . p 1 ] [ 0 ]</pre><hr>
The matrix represents the transition probabilities between m+1 states. The states are 0 in a row, 1 in a row…up to m in a row having been obtained. The columns of the matrix represent the current state, and the rows represent the next state. This is called a Markov chain. Such a solution was once presented by irchans for a similar problem.

Matrix multiplication can be carried out in Excel using MMULT to produce all of the n vectors. You have to use a ctl-shift-enter sequence to define the size of the vector as explained in the help to MMULT. An example spreadsheet is available upon request.


Derivation of exact solution 1

I give here a concise derivation of the most general form of the equation, but it is not especially easy to understand the derivation in this general form. It is easiest to refer to the derivation for the special case of m=2, p=1/2, and then insert those constants into the derivation below. It may be helpful to refer to the intial thread where the Fibonacci solution was derived for the case of m=2 and p=1/2 Fibonacci solution for coin tossing

Remember that the view we take here is that of a death and renewal process in which we consider a population of sequences at each stage, each of which “reproduces” by having a success or failure appended to it, and then some of the sequences “die” or terminate because m in a row have been obtained. We wish to determine the number that terminate at each stage which we’ll call t(n). The key realization is that this is equal to the number of sequences which do not end in a success m stages previous, which we’ll call k(n-m). This is because each of these sequences may be followed by m successes, and only those sequences will terminate at stage n. o(n) are the total number of outcomes or sequences after the nth stage. Thus:

t(n) = k(n-m)

= (1-p)*o(n-m)

= [ (1-p)/p ] * [ o(n-m-1) – t(n-m-1) ]

= [ (1-p)/p ]*{ [ 1/(1-p) ]*t(n-1) – t(n-m-1) ]

= (1/p)*t(n-1) – (1/p – 1)*t(n-m-1)

Since t(n) is the number of sequences that end in m successes at each stage, it remains to divide each of these by 2^n total sequences and sum these to get the final probability.

P(n,m,p) = sum{k = 1 to n} t(k)/2^k

In computing the above difference equation in Excel, the terms become so large that they are no longer representable. The following difference equation absorbs the division by 2^k.

F[k] = F[k-1] – (1/p – 1)*F[k-m-1]*p^(m+1) for k &gt; m+1

With initial conditions:

F[k] = 0 for k &lt; m
F[m] = p^m
F[m+1] = (1-p)*p^m

Then P(n,m,p) = sum{k = 1 to n} F[k]


Comparison of approximate and exact results

In the original thread I gave a number of approximate solutions to the problem of 3 AA in a row in 250 hands, all of varying degrees of precision and varying degrees of complexity, all of which I knew to be very useable. The lack of an exact solution led to some disparaging comments (Bruce's notion is unfortunately extremely inadequate, and from Robk Sorry Bruce, but I think you are way off on this one). I have tabulated below the results of all of the approximation methods given by me and by other posters in that thread, along with the exact answer computed by the methods described above. As can be seen, all of the approximation methods give probabilities which are quite adequate for all practical applications, and none of them are “way off” with the exception of the solution given by Robk which is indeed way off simply due to an error in his post which I have corrected. I have listed the corrected answer obtained from the method which he intended. Refer to Feller pp. 324-325. I used the goal seek feature of Excel to solve the cubic equation.

Approximations 3 and 4 appeared in the original thread. Approximations 1 and 2 are much more precise versions of 3 and 4 respectively which I have computed subsequent to that thread. The only difference is that the more precise versions use a different number of trials. Instead of 248, I have used 248 divided by the expected length of a trial which is 1*(220/221) + 2*(1/221)(220/221) + 3*(1/221)^3 = 248/1.004484216. These methods are both simple and very precise. I have given this method in other threads prior to this discussion, for example see the top of boxcars. The bottom of that post is actually an early failed attempt at a more general solution which only works for m=2. This derivation should be ignored.

The inclusion-exclusion calculation is carried out to the number of terms in the original thread. The inclusion-exclusion approximation, the Feller approximation, and the computer program all give highly precise solutions correct to the nearest whole number. The exact solutions derived above are much simpler to compute than either inclusion-exclusion or the Feller approximation.

<pre><font class="small">code:</font><hr>

Method probability (%) 1 in

1.) EXACT(both methods above) <font color="red">.002287222419 43721.15242</font color>

2.) (248/1.004484216)*(1/221)^3 .002287112016 43723.26291

3.) 1–(1-1/221)^3^(248/1.004484216) .002287319310 43719.30039

4.) 248/221^3 .002297602313 43523.63306

5.) 1 – (1-1/221)^3^248 .002297576026 43524.13103

6.) Inclusion-Exclusion .002872000000 43721.58097

7.) Cyrus’ computer program .002287224379 43721.11495

8.) Robk approximation(w/error) .002248 44483.98577

9.) Feller approximation .002287222422 43721.15236
(Robk method corrected)
1 – (1 - 1/221)*x / [ (4 – 3x)(220/221)*x^251 ]
where x is the real root of (220/221)*x*[ 1 + (1/221)*x + (1/221)^2*x^2 ] = 1
x = 1.00000092
</pre><hr>

[irchans]

Nice work Bruce. I think a general formula could be difficult to find. For the case m = 2, Mathematica gives the following formula for F

F[n] = (2^(1 - k)*p^2*(
-(1 - p - Sqrt[1 + (2 - 3*p)*p])^(-1 + k)
+(1 - p + Sqrt[1 + (2 - 3*p)*p])^(-1 + k)))/
Sqrt[1 + (2 - 3*p)*p].

For p = 1/2, that formula simplifies to

F[n] = (2^(-1 - k)*(
+ (1 + Sqrt[5]) *((1 - Sqrt[5])/2)^k
- (1 - Sqrt[5]) *((1 + Sqrt[5])/2)^k)/Sqrt[5]

= 2^(-k) * Fibonacci[k-1]

as you showed a year ago.

There is another way to approach these problems using Fourier Series. I will look that up and send this problem to a friend of mine that does combinatorics.


Cheers, Irchans

[Robk]

A few comments:

"Sorry Bruce, but I think you are way off on this one"

This was in reference to your claim that certain trials were independent (although still wrong [img]/forums/images/icons/tongue.gif[/img] ) It was not a comment on the accuracy of your (or any) approximations. Those you gave at the time were obviously excellent, given the computer simulation data which appeared in the original thread. Anyway I'm glad it motivated you.

given by Robk which is indeed way off

It's not so bad [img]/forums/images/icons/wink.gif[/img]

simply due to an error in his post which I have corrected

Thanks much for the correction. I am going to go take Feller out of the library because now I am curious about the derivation. I was working with my notes at the time which obviously sucked!

[Carl_William]

Dear Bruce:

you published or print the following on the 2+2 probability forum (for the problem in question)

” A couple of months ago, someone asked for the probability of getting 3 AA in a row in 250 hands. This gave rise to the more general question which was phrased as:

Do you have the formula for "at least 1 streak of at least m occurences in n trials"? I don't. Show us yours.: (skipping down)

Method probability (%) 1 in

1.) EXACT(both methods above) .002287222419 43721.15242

5.) 1 – (1-1/221)^3^248 .002297576026 43524.13103”



I question the answers for method 1. I feel the probablility for at least 1 streak of at least m occurences in 250 trials should be higher in method 1 than the probability in method 5. Method 5 is the binomial probability distribution for one or more occurances of a streak of three AA being dealt out for 250 sets of three consecutive holdem hands (equivalent to 750 hands). In method 5 there is no overlap between the the first set of 3 hands and the second set of three hands (and so on). Thus method 5 (in my opinion) should have a lower probability for a streak of three or more consecutive sets of two aces.

I don’t understand how to do method 1 and I don’t doubt your fine work. I just think the magnitude of the answer (probability) for method 1 is too low when compared to method 5.

Where am I wrong?

Stay well

Carl

[BruceZ]

I just think the magnitude of the answer (probability) for method 1 is too low when compared to method 5.
Where am I wrong?

First of all, there is a typo in method 5, and in method 3 as well. The numbers are correct, but the expressions should be:

method 5: 1 - (1 - 1/221^3)^248
method 3: 1 - (1 - 1/221^3)^(248/1.004484216)

The way I had it would give the probability of getting a single AA in 3 times 248 hands which would be very high. You probably still think that the exact answer of method 1 should be higher than method 5, but it is correct that it is lower.

Method 5 computes 1 minus the probability of no runs of length m or greater. The probability it computes of no runs is smaller than the exact number, so the final answer will be greater than the exact answer. The reason method 5 computes the probability of no runs smaller than exact is because it treats the probabilities of no run starting on each hand as independent and equal to 1 - 1/221^3 per hand. The exact result would be:

P(no run) = P(1)*P(2|1)*P(3|1,2)*...P(248|1,2,...247)

Where P(k) = P(run NOT starting on k) and P(k|j) = P(run NOT starting on k | run NOT starting on j).

This (exact) probability will be larger than 1 - (1/221^3)^248. For example, P(1) = 1 - 1/221^3. P(1|2) &gt; 1 - 1/221^3. The reason is that if a run DID start on hand 1, the probability of a run starting on hand 2 would be just the probability of AA on hand 4 since we already know we have AA on hands 2 and 3, and this probabity is just 1/221 &gt; 1/221^3. Thus the probability of a run starting on hand 2 given that one did NOT start on hand 1 or must be &lt; 1/221^3, and the probability of a run NOT starting on hand 2 given that a run did NOT start on hand 1 or P(1|2) would be &lt; 1/221^3 or &gt; 1 - 1/221^3. This will be true for all the terms, so the overall product in method 5 is smaller than the exact method, so it produces a probability of a run which is larger than the exact method.

I also originally thought that method 5 should produce a lower bound, but I then realized even before I derived the exact answer that it would in fact be an upper bound, and I posted this here .

Method 1 must give the exact answer because it gives the identical answer as the Markov matrix method which must be correct. I also ran both my formula and the Markov matrix for m=5, n=250, p=1/221 and the results were again identical.

[BruceZ]

This was in reference to your claim that certain trials were independent (although still wrong)

That's correct, the trials are not independent; however, to a high degree of approximation we can treat them as being independent. If this were not the case, approximations 5 and even better 3 would not be close by definition. I had reworded this claim in the original thread, but it did not appear until after you had posted.

I am going to go take Feller out of the library because now I am curious about the derivation. I was working with my notes at the time which obviously sucked!

You just added q in the denominator where it should have been multiplied.

Feller derives a difference equation in a different variable than I did, but in a footnote he alludes to the fact that the classical approach derives one in the same variable that I used. I didn't notice this note until after I had done my derivation. Feller never mentions the Fibonacci solution even though he has data on coin tossings. It is possible that he was not aware of this result. From his remarks it also seems possible he had not seen the general form that I have derived here.

Feller's solution is approximate as it only uses one root of the polynomial. This root is dominant, so the approximation is a good one. To make his method exact, all of the roots must be found, and this is much more complicated. So while his footnote claims that my approach is more complicated, it is much simpler for obtaining an exact solution.

[BruceZ]

irchans,

I think a general formula could be difficult to find. There is another way to approach these problems using Fourier Series.

My solution is general in the sense that it gives a solution for all m,n,p and all the necessary terms are straightforward to compute. By general you mean F[n] = with no more reference to F[n] (no recursion). For that you would need to take the z-tranform (or generating function) of the difference equation, solve for F(z) = 1/(mth order polynomial in z), find all the roots of the polynomial, do a partial fraction expansion, and then invert the individual terms back to the n domain to get a sum of exponentials in n. For m=2, p=1/2 the roots are [ 1 +/- sqrt(5)] /2 as you found. It turns out that the real root always dominates, and this is what Feller's solution uses.

Is this what you mean by the Fourier series method? The z-transform is a general case of the discreet Fourier tranform. This is the technique I'm very familiar with.

I take it that mathematica doesn't spit out a result for general m? No closed form solution exists for factoring a general mth order polynomial, but it is not impossible that this particular family of polynomials could have a general solution.

I have found a non-recursive solution for F using only the single real root and a single term of F[n]. It is analogous to the reduced form of the Fibonacci solution.

BTW, since you will be looking at this closely, this line of the derivation:

= [ (1-p)/p ]*{ [ 1/(1-p) ]*t(n-1) – t(n-m-1) ]

uses a substitution from previous lines. This wasn't too obvious.

-Bruce

[Carl_William]

Please disregard my post on July 15,2003 concerning "Re: Great problem solved - original material." I have thought things over and have a new slant which I think is pertinent. I will shortly post my new slant.

thank you

[Carl_William]

Re: Great problem solved - original material – My new slant….

Dear Bruce and others interested in these types of problems,

In elementary probability problems where the probabilities are constant (Bernoulli trials), the use of the Markoff process is not necessary. Like they say: “in fixed or constant probability problems” – the cards, dice, and wheel, “whatever” has no memory. Therefore in a set of three deals of holdem hands, the odds are fixed and never change. For the problem in question: the answer can be determined by the elementary Bernoulli trial formula of sampling with replacement:

P = comb(n,x) * p^x * (1-p)^(n-x)

Where n = 248, p =(1/221)^3, x = number of successes,
and as you posted in Method 5: Answer = 1 – (1-1/221)^3^248

This model used here can be thought of as: 248 sets of three deals (or 3*248 = 744 deals) – each set taken as independent of any other set. The probability is constant so this is a valid assumption. Look at it this way: set 1 = deal 1 + deal 2 + deal 3 – this is one trial; set 2 = deal 4 + deal 5 + deal 6; …. ; last set 248 = deal 742 + deal 743 + deal 744. Like I am trying to say, “this is valid because the cards have no memory,” and simplifies to the method 5 formula: Answer = 1 – (1-1/221)^3^248.

If Methods(1) are also exact, then they should calculate the same answer as method 5. As you know, the choice of 250 deals: 248 sets of 3 deals is purely arbitrary – also the selection of p(AA) = 1/221 is also arbitrary:for Answer = 1 – (1-1/221)^3^248. A simple example problem with flipping a coin where two consecutive heads appears in say 4 trials could easily verify that the simple Bernoulli trial formula works…. This simple problem can be enumerated on an Excel spreadsheet. I would check the boundary (end conditions for method(s) 1 – maybe they are tricky to handle.

Markoff process are generally use in dynamic programming problems where the odds change with time, and recent and past history can be used to determine future actions.

This is just my opinion and hopefully I don’t have too many mistakes….

Regards &amp; stay well

Carl

[Carl_William]

Excuse my typo error:
I repeated or carried on the error – should be:

Method 5.) Answer = 1 – (1-1/221^3)^248 =2.2975760259003E-05

not: 1 – (1-1/221)^3^248

---------------
Re: Great problem solved - original material – My new slant…. I think I understand the problem definition.

Dear Bruce and others interested in these types of problems,

In elementary probability problems where the probabilities are constant (Bernoulli trials), the use of the Markoff process is not necessary. Like they say: “in fixed or constant probability problems” – the cards, dice, and wheel, or “whatever” have no memory. Therefore in a set of three deals of holdem hands, the odds are fixed and never change. For the problem in question: the answer can be determined by the elementary Bernoulli trial formula of sampling with replacement:

P = comb(n,x) * p^x * (1-p)^(n-x)

Where n = 248, p =(1/221)^3, x = number of successes,
and as you posted in Method 5: Answer = 1 – (1-1/221^3)^248
=2.2975760259003E-05

MY LOGIC IS:
This model used here can be thought of as: 248 sets of three deals (or 3*248 = 744 deals) – each set taken as independent of any other set. The probability is constant so this is a valid assumption. Look at it this way: set 1 = deal 1 + deal 2 + deal 3 – this is one trial; set 2 = deal 4 + deal 5 + deal 6; …. ; last set 248 = deal 742 + deal 743 + deal 744. Like I am trying to say, “this is valid because the cards have no memory,” and simplifies to the method 5 formula: Answer = 1 – (1-1/221^3)^248.

If Methods(1) are also exact, then they should calculate the same answer as method 5. As you know, the choice of 250 deals: 248 sets of 3 deals is purely arbitrary – also the selection of p (AA) = 1/221 is also arbitrary: for Answer = 1 – (1-1/221^3)^248. A simple example problem with flipping a coin where two consecutive heads appears in say 4 trials could easily verify that the simple Bernoulli trial formula works…. This simple problem can be enumerated on an Excel spreadsheet. I would check the boundary (end) conditions for method(s) 1 – maybe they are tricky to handle.

Markoff processes are generally used in dynamic programming problems where the odds change with time, and recent and past history can be used to determine future actions. Also, in the real world, many answers that are within 1% or 5% are adequate -- no always but usually.


These are just my opinion and hopefully I don’t have too many mistakes.

Regards &amp; stay well

Carl