B Probability that 1000 coin flips results in >600 tails

[Jonathan212]

What is the formula for this? Throwing a flipping coin N times, what is the probability that the number of tails results is higher than M?

 

[BvU]

Well, doesn't friend google answer your question ? What did you find ?

 

[Jonathan212]

This is not homework if that's what you were thinking. Sounds specialized enough to me that a google answer is not readily available. Unless you dig a lot in which case might as well derive it from scratch.

 

[PeroK]

This is not homework if that's what you were thinking. Sounds specialized enough to me that a google answer is not readily available. Unless you dig a lot in which case might as well derive it from scratch.

Try "Binomial distribution" and the "Normal approximation to the binomial distribution".

 

[BvU]

"flip coin probability"does the trick. That's not 'digging in deep', really

 

[Jonathan212]

Thanks. So the answer is:

Pr(M; N, 0.5) = N! / (M! * (N - M)!) * 0.5^M * (1 - 0.5)^(N-M)

https://en.wikipedia.org/wiki/Binomial_distribution

Is this High School level in the US nowadays?

 

[PeroK]

Thanks. So the answer is:

Pr(M; N, 0.5) = N! / (M! * (N - M)!) * 0.5^M * (1 - 0.5)^(N-M)

https://en.wikipedia.org/wiki/Binomial_distribution

Is this High School level in the US nowadays?

Note that $(1 - 0.5) = 0.5$, so you can simplify that expression. That is for precisely $M$ tails in $N$ flips. To get the answer to your original question you have to take the sum from $M = 601$ to $M = 1000$.

Note also that this is where the normal approximation is useful. Then, you can get the answer by looking up a normal distribution table.

 

[Jonathan212]

Oopsa. Excel can't deal with N > 170. Any trick to avoid the overflow of "171!" ?

 

[PeroK]

Oopsa. Excel can't deal with N > 170. Any trick to avoid the overflow of "171!" ?

Normal approximation!

https://revisionmaths.com/advanced-level-maths-revision/statistics/normal-approximations

 

[Jonathan212]

Got stuck now. Something is wrong. Do we want a normal distribution with a standard deviation of N * 0.5 * ( 1 - 0.5 ) and a mean of N * 0.5? Excel has the NORMDIST() function but it returns a result 4 times bigger than the binomial distribution at N = 170, M = 102.

 

[PeroK]

Got stuck now. Something is wrong. Do we want a normal distribution with a standard deviation of N * 0.5 * ( 1 - 0.5 ) and a mean of N * 0.5? Excel has the NORMDIST() function but it returns a result 4 times bigger than the binomial distribution at N = 170, M = 102.

https://www.statisticshowto.datasci...theorem/normal-approximation-to-the-binomial/

 

[Jonathan212]

So we want sigma = sqrt( N * 0.5 * ( 1 - 0.5 ) ). And the answer for N = 1000, M = 601 is 1 in 21 million. Hurray! That's statistically significant as hell.

Looking at the error of normal versus the binomial up to N = 170 (M = floor(0.6*N+0.5)):

2015 sti 0 60

It actually gets worse after N = 140. Is it meant to?

 

[PeroK]

So we want sigma = sqrt( N * 0.5 * ( 1 - 0.5 ) ). And the answer for N = 1000, M = 601 is 1 in 21 million. Hurray! That's statistically significant as hell.

Looking at the error of normal versus the binomial up to N = 170:

View attachment 238253
2015 sti 0 60

It actually gets worse after N = 140. Is it meant to?

I'm taking a guess that the normal distribution, being essentially $e^{-x^2}$, will fall off to zero faster than the binomial. So, after a certain point the relative error will increase, but relative to some very small probabilities.

 

[Jonathan212]

At N = 100, M = 51, I get 0.078 from both the binomial and the normal. Sure this is right?

 

[PeroK]

At N = 100, M = 51, I get 0.078 from both the binomial and the normal. Sure this is right?

For N =100, I would expect almost all cases to lie in the interval 40-60 tails. That's about 5% each on average. So, 8% chance of exactly 51 tails sounds about right. Don't you think?

 

[Jonathan212]

But I want M or more tails, not M.

 

[Jonathan212]

Wikipedia: "the binomial distribution with parameters n and p is the discrete probability distribution of the number of successes in a sequence of n independent experiments".

 

[PeroK]

But I want M or more tails, not M.

Then you have to add all cases from M = 51 to 100. See post #7.

Your previous answer in post #12 may be for exactly 601 tails.

 

[Jonathan212]

Looks like there is no formula that saves you having to do the sum of the binomial. But luckily, there is for the normal in excel.

 

[PeroK]

Looks like there is no formula that saves you having to do the sum of the binomial. But luckily, there is for the normal in excel.

There should be a "cumulative" field for both.

 

[BvU]

You want to be careful with this 'cumulative' output: BINOM.DIST(50,100,0.5,TRUE) gives you 0.54 , not 0.5 !

 

[Jonathan212]

Back to basics cause I don't really know what I'm doing.

N = 3 throws, possible outcomes: 000, 001, 010, 011, 100, 101, 110, 111

0 heads -> p = 1/8
1 heads -> p = 3/8
2 heads -> p = 3/8
3 heads -> p = 1/8

Binomial distribution in Excel is as expected (NO cumulative flag):

0 heads -> p = 0.125
1 heads -> p = 0.375
2 heads -> p = 0.375
3 heads -> p = 0.125

Binomial distribution in Excel with cumulative flag:

up to 0 heads -> p = 0.125
up to 1 heads -> p = 0.5
up to 2 heads -> p = 0.875
up to 3 heads -> p = 1

We want at least M heads:

M = 0 -> p = 1
M = 1 -> p = 7/8
M = 2 -> p = 4/8
M = 3 -> p = 1/8

Now how do we write the last list in terms of the binomial function of Excel? As follows it is WRONG:

= 1 - BINOMDIST( M , N , 0.5, 1 )

 

[BvU]

= 1 - BINOMDIST( M-1 , N , 0.5, 1 )

[edit] unfortunately, the function isn't defined for M-1 = -1

 

[mfb]

=IF(M=0,1,1 - BINOMDIST( M-1 , N , 0.5, 1 ))

Ugly but it works.

 

[StoneTemplePython]

I'm taking a guess that the normal distribution, being essentially $e^{-x^2}$, will fall off to zero faster than the binomial. So, after a certain point the relative error will increase, but relative to some very small probabilities.

qualitatively this can't be the whole story -- the binomial is bounded while the Gaussian is unbounded...

 

[StoneTemplePython]

@Jonathan212

what are you ultimately looking for here? Something to type into excel? If you're doing it computationally, you maybe be able to just sum over 400 terms or so. If you know basic coding, there's a slick way to do this with trees. If you want something 'nice' as an approximation, that seems ok. Normal/ Gaussian approximation works though there aren't easy closed form integrals (which is what you'd want).

You can get easy upper and lower bounds on Gaussian integrals though. It may also be worthwhile to know that asymptotically for a standard normal (i.e. mean zero, std deviation of 1: so you'd need to scaled and shift your problem to 'become standard' to use this) you have, for some $c \gt 0$
$P(X\gt c) \sim \frac{1}{c\sqrt{2 \pi}} e^{-\frac{c^2}{2}}$
(which is a very nice approximation for $c \gt 1$ -- technically it is an upper bound on the integral which has a simple proof using the 'slip-in' trick)

the bound is below

Spoiler

for $c \gt 0$ and standard normal r.v. $X$, we have
$P(X\gt c) = \frac{1}{\sqrt{2 \pi}}\int_c^\infty e^\frac{-x^2}{2} dx$
$= \frac{1}{\sqrt{2 \pi}}\int_c^\infty \frac{x}{x} e^\frac{-x^2}{2} dx$
$\leq \frac{1}{\sqrt{2 \pi}}\int_c^\infty \frac{x}{c} e^\frac{-x^2}{2} dx$
$= \frac{1}{c\sqrt{2 \pi}}\int_c^\infty x e^\frac{-x^2}{2} dx$
$= \frac{1}{c\sqrt{2 \pi}} e^{-\frac{c^2}{2}}$

learning how to convert a binomial into a normal and then how to convert a normal into a standard normal, seems like a worthwhile goal for this
- - - -
Also, there are relatively straightforward proofs of the Chernoff Bound for coin tossing...

 

[Jonathan212]

what are you ultimately looking for here?

The way to assess the statistical significance of a true random number generator's observed bias. Got the answer in Excel, it is as BvU says. Only thing is, the result is a bit suspicious for high N: probability of 600 OR MORE heads in 1000 throws is 1 in 7.3 billion.

 

[PeroK]

The way to assess the statistical significance of a true random number generator's observed bias. Got the answer in Excel, it is as BvU says. Only thing is, the result is a bit suspicious for high N: probability of 600 OR MORE heads in 1000 throws is 1 in 19 billion.

Getting 600 or more heads in 1000 flips sounds next to impossible to me. The probablity of getting 60 or more in 100 flips is about 3%. And, more or less, you'd have to do that 10 times in a row.

 

[Jonathan212]

Look at some more numbers. All at 60% heads or more:

Probability of 15 or more heads in 25 throws = 1 in 4.7
Probability of 60 or more heads in 100 throws = 1 in 35
Probability of 150 or more heads in 250 throws = 1 in 1,062
Probability of 240 or more heads in 400 throws = 1 in 27,000
Probability of 300 or more heads in 500 throws = 1 in 220,000
Probability of 360 or more heads in 600 throws = 1 in 1,800,000
Probability of 480 or more heads in 800 throws = 1 in 1,200,000,000

 

[PeroK]

Look at some more numbers. All at 60% heads or more:

Probability of 15 or more heads in 25 throws = 1 in 4.7
Probability of 60 or more heads in 100 throws = 1 in 35
Probability of 150 or more heads in 250 throws = 1 in 1,062
Probability of 240 or more heads in 400 throws = 1 in 27,000
Probability of 300 or more heads in 500 throws = 1 in 220,000
Probability of 360 or more heads in 600 throws = 1 in 1,800,000
Probability of 480 or more heads in 800 throws = 1 in 1,200,000,000

That's the law of averages for you!

 

[Jonathan212]

Alright but such a combinatorial explosion was not expected.

 

[StoneTemplePython]

The way to assess the statistical significance of a true random number generator's observed bias. Got the answer in Excel

'true random number generator' is a bit of a loaded term. Excel isn't known to be good for pseudo random number generators.

Got the answer in Excel, it is as BvU says. Only thing is, the result is a bit suspicious for high N: probability of 600 OR MORE heads in 1000 throws is 1 in 19 billion.

you're could run into numeric stability issues if you're calculating and summing extremely small probabilities so extra care is needed. It sounds about right assuming a fair coin. I know the correct probability is cannot be bigger than $\frac{2.07}{10^{9}}$. Here's why:

you know every probability $p$ satisfies $0\leq p \leq 1$. Now using chernoff bounds, I can say, with

$S_n = X_1 + X_2 + ... + X_n$
where each $X_k$ is a indendent bernouli (coin toss) with value 1 for tails and 0 for heads.
i.e. $S_n$ is a random variable that counts the number of tails

So the expected number of tails $= \mu = E[S_n] = n \cdot E[X_k] = 1000 \cdot 0.5 = 500$

what the Chernoff bound tells us is

$0 \lt Pr\Big\{\text{600 or More Tails}\Big\} = Pr\Big\{\big(S_n - \mu\big) \geq 100\Big\} \leq \exp\big(\frac{-2\cdot 100^2}{1000}\big) \lt 2.07 \cdot 10^{-9}$

(note thread title says $\gt 600$ tails where as a more recent posting says $\geq 600$, I went with the latter)
- - - - -

edit:
for a basic look at Chernoff Bounds for coin tossing, consider looking at pages 20-30, here:
https://ocw.mit.edu/courses/electri...ce-fall-2010/readings/MIT6_042JF10_chap19.pdf

part of the math for CS course, here:
https://ocw.mit.edu/courses/electri...tics-for-computer-science-fall-2010/readings/

 

[Jonathan212]

Indeed Excel's cumulative binomial runs out of steam at N = 1029, it fails above that value. But so does Excel's normal distribution at N = 1555 (the cumulative option). Btw the source of my true random data is not Excel but /dev/random or "haveged".

 

[Jonathan212]

The plot looks alright though.

X-axis = N
Y-axis = 1 / probability that N throws result in >=N*0.6 heads.

 

[Matt Benesi]

It's 100% if you do 1000 coin flips again and again until you get greater than 600 tails.

Well, unless you're cursed. Then it's 0%, even with a double tailed coin. You get killed before you flip the coin the 1000th time.

 

[StoneTemplePython]

Look at some more numbers. All at 60% heads or more:

Probability of 15 or more heads in 25 throws = 1 in 4.7
Probability of 60 or more heads in 100 throws = 1 in 35
Probability of 150 or more heads in 250 throws = 1 in 1,062
Probability of 240 or more heads in 400 throws = 1 in 27,000
Probability of 300 or more heads in 500 throws = 1 in 220,000
Probability of 360 or more heads in 600 throws = 1 in 1,800,000
Probability of 480 or more heads in 800 throws = 1 in 1,200,000,000

In some fundamental sense, the problem is that you are looking at this the wrong way -- you are looking at this in terms of some kind of constant percentage -- i.e. 60% of n.

A much improved way to look at this is consider your threshold (e.g. at least 15 heads ) vs the 'distance' from that to the mean. The distance involved is called the amount of standard deviations from the mean. If you re-visit these numbers, you'll see that you're progressively increasing the distance between the threshold and the mean as you go down this table. Look at a picture of the bell curve and you can see how it rapidly drops off the farther from the mean you go... the result here should not surprise you if you think about it this way (and know/assume that a normal approximation is reasonable).

So, here's quick look:

$\left[\begin{matrix} \text{n} & \geq\text{threshold} & \text{raw excess over mean} & \text{ std deviations from the mean} &\text{normal approximation} \\ 25 & 15 & 2.5 & 1 & 0.24 \\ 100 & 60 & 10 & 2 & 0.027 \\ 250 & 150 & 25 & 3.16 & 0.00085 \\ \vdots& \vdots& \vdots & \vdots & \vdots \\ 500 & 300 & 50 & 4.47 & 4.0 \cdot 10^{-6} \\ 1000 & 600 & 100 & 6.32 & 1.3 \cdot 10^{-10} \\\end{matrix}\right]$

for a binomial, you have
$\text{mean}= \mu = n \cdot p$ and
$\text{variance}= n \cdot p(1-p) \to \text{standard deviation}= \sqrt{n \cdot p(1-p)}$

where the normal approximation used was:
$P(X\geq c) = P(X\gt c) \sim \frac{1}{c\sqrt{2 \pi}} e^{-\frac{c^2}{2}}$
(reference post #26 for more information)

and $c$ is precisely the number of standard deviations from the mean that your threshold is.

(Note: we can fine tune this with the 1/2 / continuity correction as is common for approximating binomials with standard normals, among other techniques though I skipped it them for convenience. A more meticulous approach would also bound the error made when approximating a binomial by a normal. )

 

[Jonathan212]

Does the integral of the standard distribution have an analytical solution by any chance? If not, how do you know what interval to use to the 2 significant figures in the probability

 

[BvU]

It has a name: 'error function' and there is an expression for it, but no analytical solution.

There are plenty detailed tables.

Can't say I understand your second question, though.

 

[Jonathan212]

Second question was about generating the look up table by numerical integration. If 2 s.f. are wanted for the probabilities, 2 s.f or more are required for the look up table, but an arbitrary integration step of, say, sigma*0.1 does not guarantee 2 s.f. or does it?

 

[StoneTemplePython]

Does the integral of the standard distribution have an analytical solution by any chance? If not, how do you know what interval to use to the 2 significant figures in the probability

it's a B thread so I won't push this too hard: once you get away from really simple problems you won't get nice analytical solutions. Better to approximate and bound.

I can tell you that for a standard normal $X$ and any $c\gt 0$
$\big(\frac{1}{{\sqrt{2\pi}}}e^{-\frac{c^2}{2}}\big)\big(\frac{c}{c^2 + 2}\big) \leq P(X\geq c) \leq \big(\frac{1}{\sqrt{2 \pi}} e^{-\frac{c^2}{2}}\big) \big(\frac{1}{c}\big)$

note that the ratio of the upper and lower bound tend to one. Thus you can say that asymptotically
$P(X\geq c) \sim \big(\frac{1}{\sqrt{2 \pi}} e^{-\frac{c^2}{2}}\big) \big(\frac{1}{c}\big)$

which is what I used in that table in post #36.

I proved the upper bound in post #26. Unfortunately all of the various asymptotically tight lower bounds that I'm aware of are a bit more involved.

edit:
here's a classical way of getting a pretty good bound, courtesy of Feller

where $\phi(x) := \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$
(i.e. the density for the standard normal random variable)

consider that for any $x \gt 0$ we have the point wise bound

$\big(1- 3x^{-4}\big)\phi(x) \lt \phi(x) \lt \big(1+x^{-2}\big)\phi(x)$

now integrating over this bound, for any $c \gt 0$
$\big(c^{-1} - c^{-3}\big)\phi(x) = \int_{c}^\infty \big(1- 3x^{-4}\big)\phi(x)dx \lt \int_{c}^\infty\phi(x)dx = P\big(X\gt c\big) \lt \int_{c}^\infty\big(1+x^{-2}\big)\phi(x)dx = \big(c^{-1}\big)\phi(x)$

or if you prefer:
$\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\big(\frac{1}{c} -\frac{1}{c^3}\big) \lt P\big(X\gt c\big) \lt \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\big(\frac{1}{c}\big)$

where the upper bound is the same as before

- - - -
edit:
if you want to get the numerical values of the Gaussian integral, can't you just use a built-in excel function?

Indeed Excel's cumulative binomial runs out of steam at N = 1029, it fails above that value. But so does Excel's normal distribution at N = 1555 (the cumulative option).

or maybe not... though this does seem a bit odd to me. There really should not be an "N" input for the normal distribution. I think you are doing something wrong.

also

Btw the source of my true random data is not Excel but /dev/random or "haveged".

whatever this means, it sounds like you know how to program. Certainly I've looked up CDF (or complementary CDF in this case) values for standard normal random variables using Scipy on many occasions, which would be my suggestion. I guess I'm scratching my head a bit as to why you're using Excel, but I digress.

for example:

from scipy.stats import norm
print(1-norm.cdf(4.47))
print(1-norm.cdf(6.32))

prints out "3.9109798603e-06" and "1.30781607766e-10" which are basically the same thing as the bottom right two entries of the table post #36. Notice that there are no "N" values here.

This result should either give you confidence in scipy or the asymptotic approximations of the Gaussian integral that I used (or maybe increase confidence in both).

 

[Jonathan212]

if you want to get the numerical values of the Gaussian integral, can't you just use a built-in excel function?

Indeed. Looks like the thing to do is generate the look up table with Excel for maximum accuracy, save it as a file, and load the file from a C or PHP program that handles the TRNG source to simulate lots of coin flips per second and look up the resulting count of tails in the table, in order to display the statistical significance of the result every second. Guess what comes next.

Need some help with Excel's function. Got the following formula for the integral but I think it is wrong, it does not match the binomial sum at all.

= 1 - NORMDIST( M - 1, N * 0.5, SQRT( N * 0.5 * (1-0.5) ), 1 )

Correct binomial sum is:

= 1 - BINOMDIST( M - 1, N, 0.5, 1 )

These are the formulas where Excel fails from N = 1555 upwards for the normal, and from N = 1030 upwards for the binomial.

EDIT: could also write the normal as follows and still get the same numbers while it now fails at N = 1718:

= NORMDIST( N - M, N * 0.5, SQRT( N * 0.5 * (1-0.5) ), 1 )

At N = 1000 it still has a 6.5% error relative to the binomial.

 

[StoneTemplePython]

At N = 1000 it still has a 6.5% error relative to the binomial.

To be clear I believe you mean $\frac{\big \vert\text{binomial probability} - \text{normal approx}\big \vert}{\text{binomial probability}} = 6.5\%$

This too should not surprise you. For the N=1000 case you are talking about a vanishingly small probability of having your event of $\geq 600$ tails occur (which has a probability less than $\frac{2.07}{10^{9}}$), so even if the normal is almost a perfect match to the binomial over the entire distribution (total variation distance is the underlying concept I'm hinting at) any very small differences in calculated values will be magnified when talking about extremely rare events. Put differently, it is well understood that normal approximations are weakest when you go deep into the tails / rare events. Having a $\approx 94\%$ match for something this rare is surprisingly good.

This is also part of the reason I suggested using Chernoff Bounds in post 32 -- you can get an exact upper bound and not have to worry about approximation errors for rare events.

 

[slappmunkey]

The probability that you flip a coin 1000 times and 600 times is tail?
50%...
That's the same question just reworded as "if I flip a coin 10 times and 9 times it's heads, what's the odds it will land tails on the 10th flip?" , Same answer, 50%.

 

[mfb]

The question is about seeing (at least) 600 tails in total, not about the probability of a specific flip out of the 1000.

That's the same question just reworded as "if I flip a coin 10 times and 9 times it's heads, what's the odds it will land tails on the 10th flip?"

It isn't.

 

[slappmunkey]

The question is about seeing (at least) 600 tails in total, not about the probability of a specific flip out of the 1000.It isn't.

It is. Each individual flip is a 50% chance to land either way, that means it's a 50/50% chance no matter how many flips you do.

 

[StoneTemplePython]

It is. Each individual flip is a 50% chance to land either way, that means it's a 50/50% chance no matter how many flips you do.

two ideas:
i) graphical: pick any row n ($n\geq 2$) of pascal's triangle. You can use row 1000 if you want. Are the numbers all the same? This has something to do with the binomial distribution.
ii) suppose for a contradiction that what you've said is true and we toss a coin 1000 times. Then the probability of 0 tails is $\frac{1}{2}$, probability of 1 tails is $\frac{1}{2}$, probability of 3 tails is $\frac{1}{2}$, ..., probability of 999 tails is $\frac{1}{2}$, and probability of 1000 tails is $\frac{1}{2}$. So the total probability is $\frac{500500}{2}\gt 1$ which is contradiction because probabilities sum to one.

 

[slappmunkey]

You are over thinking it.

Break it down flip by flip.
Flip 1: 50%
Flip 2: 50%
Flip... 1000: 50%

Now what's the average of all the flips... 50%...

 

[fresh_42]

The OP's question has been answered.

Thread closed.