# Emperically testing a "p-test" (of multiple fair coin flips)

In a previous question, I asked how to test if a coin is fair. Now I want to empirically test if this test works.

One answer is that programs like R and python have built-in "binomial p-tests" that can be called to do this.

Here is an example of some python code to do such a p-test for a single case of flipping 1000 fair coins:

import numpy as np
from numpy import random
import scipy
from scipy import stats

def flipAndDoPTest(numberOfFlips, weight):
flippedCoinSimulation = np.random.binomial(1, weight, numberOfFlips) #first input changes sum of coins
numberOfTails = np.sum(flippedCoinSimulation==0)
return pvalue

numberOfFlips = 1000
weight = .50
ptestvalue = flipAndDoPTest(numberOfFlips, weight)
if ptestvalue>.05:
print("the ptest has a value of:", ptestvalue)
print("The null hypothesis cannot be rejected at the 5% level of significance because the returned p-value is greater than the critical value of 5%.")
if ptestvalue<.05:
print("the ptest has a value of:", ptestvalue)
print("The null hypothesis can be rejected at the 5% level of significance because the returned p-value is less than the critical value of 5%.")


Now I want to empirically test what this "5% level of significance" means. It seems like there is a lot of disagreement over the interpretation of p-values, so I want to just simulate what occurs in my case.

First I would like to test if 5% of the time a p-value less than .05 will occur for a fair coin. To do this, I repeat this p-test 1000 times (and each p-test is for the event of flipping a fair coin 10000 times). Now I collect all of the times the p-value is less than .05. The code for this is here:

numberOfFlips = 10000
weight = .50
numberOfTests = 1000
StatisticalSignificanceProbability = .05
pTestList = np.zeros(numberOfTests) #initialization
for i in range(numberOfTests):
#for each i in the loop, do a p-test of 10,000 fair coin flips and add it to a list
ithPTest = flipAndDoPTest(numberOfFlips, weight)
pTestList[i] = ithPTest
#take this list and count all of the times there are cases below .05
numberOfSheerCoincidences = sum(pTestList<StatisticalSignificanceProbability)
expectedNumberOfSheerCoincidences = numberOfTests*StatisticalSignificanceProbability

print("numberOfSheerCoincidences: ", numberOfSheerCoincidences)
print("expectedNumberOfSheerCoincidences: ", expectedNumberOfSheerCoincidences)


Now I expect that 5% of my 1000 p-tests are going to be less than .05 (so .05*1000 = 50). But every time I run it I get a number that is significantly smaller than 50. Now this result has a random distribution, so I then wrote code to repeat this procedure to get a histogram distribution of the resultant data:

numberOfFlips = 100
weight = .50
numberOfDataPoints = 1000
pTestResultsDataPoints = np.zeros(numberOfDataPoints) #initialization
for j in range(numberOfDataPoints):
#repeating this collection of p-test to get a range of different values
numberOfTests = 1000
StatisticalSignificanceProbability = .05
pTestList = np.zeros(numberOfTests) #initialization
for i in range(numberOfTests):
ithPTest = flipAndDoPTest(numberOfFlips, weight)
pTestList[i] = ithPTest
numberOfSheerCoincidences = sum(pTestList<StatisticalSignificanceProbability)
expectedNumberOfSheerCoincidences = numberOfTests*StatisticalSignificanceProbability
pTestResultsDataPoints[j] = numberOfSheerCoincidences

n, bins, patches = plt.hist(pTestResultsDataPoints, 50)
plt.show()


And with that result I get a distribution that is centered around 35 and not 50.

Is this result expected? I was expecting a normal distribution around 50.

In general, there is no such binomial test with significance level exactly $$\alpha = 0.05,$$ on account of the discreteness of binomial distributions.

For an exact test at $$\alpha = 0.05$$ based on a continuous test statistic, the distribution of P-value when $$H_0$$ is true would be standard uniform and the probability that the P-value is below $$0.05$$ would be exactly $$0.05.$$

If $$n = 100,$$ testing $$H_0: p = .5$$ against $$H_a: p \ne 0.5,$$ the closest one can get to a test at the 5% level (without going over 5%) is $$0.0352 = 3.52\%.$$

2*(1 - pbinom(60, 100, .5))
[1] 0.0352002
2*(1 - pbinom(59, 100, .5))
[1] 0.05688793


[It does not help to use a normal approximation with nominal 5% level because z values close to $$\pm 1.96$$ cannot be achieved. For $$\mathsf{Binom}(100,.5)$$ the normal approximation is very accurate so it hardly matters whether one does an 'exact' binomial test or an approximate normal one.]

Below I simulate 100,000 tests in R for $$n=100$$ observations, and summarize and plot histograms of P-values. Probabilities should be accurate to about two places. The binomial test 'binom.test' and approximate normal test 'prop.test' both have the anticipated P-values.

set.seed(2021);  n = 100
pv.b = replicate(10^5,
binom.test(rbinom(1,n,.5),n,.5)$p.val) mean(pv.b < 0.05) [1] 0.03605 # aprx 0.0352 2*sd(pv.b < 0.05)/sqrt(10^5) [1] 0.001178995 # aprx 95% margin of sim error set.seed(2021); n = 100 pv.n = replicate(10^5, prop.test(rbinom(1,n,.5),n,.5)$p.val)
mean(pv.n < 0.05)
[1] 0.03605


The figures below show the simulated distributions of the P-values under $$H_0.$$

R code for figure:

par(mfrow=c(1,2))
hist(pv.b, prob=T, xlim=c(-.01,1.01), col="skyblue2")
abline(v = .05, col="red")
hist(pv.n, prob=T, xlim=c(-.01,1.01), col="skyblue2")
abline(v = .05, col="red")
par(mfrow=c(1,1))


Notes: (1) At the resolution of the plot, the two histograms look the same, but there are miniscule differences in the distributions; at the borderline the binomial and approximate normal tests had slightly different P-values in 8070 tests out of 100,000. But they never disagreed about rejection at the 5% level.

sum(pv.b == pv.n)
[1] 8070
mean(abs(pv.b - pv.n))
[1] 0.0001792384
sum((pv.n <= .05) == (pv.n <=.05))
[1] 100000


(2) Without using a randomized test, the closest $$\alpha$$ to 5% without going over for $$n=1000$$ is $$0.046 = 4.6\%,$$

2*(1-pbinom(531,1000,.5))
[1] 0.0462912
2*(1-pbinom(530,1000,.5))
[1] 0.05367785


It is (nearly) impossible for a test statistic that has a discrete distribution to attain exactly a 5% (or any other) significance level. So instead, many so-called "exact" tests use the largest attainable significance level that is less than $$\alpha$$ (or 5% in your case). That is what you are seeing.

Here is a reference that elaborates:

https://www.jstor.org/stable/2684683?origin=crossref