resampling.py

from math import factorial as fact
from random import random, randrange, shuffle
from statistics import mean

def combinations(n, r):
    return fact(n) / fact(r) / fact(n-r)

def cum_binom(n, r, ph):
    total = 0.0
    for h in range(r, n+1):
        total += ph ** h * (1.0 - ph) ** (n - h) * combinations(n, h)
    return total

# Estimate the probability of getting 42500 or more heads
# from 70000 spins of biased coin (60% heads and 40% tails)    

# Estimate the probability of getting 5 or more heads
# from 7 spins of biased coin (60% heads and 40% tails)
print(cum_binom(7, 5, 0.60))

def trial():
    nh = 0
    for i in range(7):
        nh += random() < 0.60
    return nh >= 5

cases = 0
for i in range(10000):
    cases += trial()
print(cases / 10000.0)

trial = lambda : sum(random() < 0.60 for i in range(7)) >= 5
print(sum(trial() for i in range(10000)) / 10000.0)

# Probability of the median of 5 samples being in the
# middle two quartiles of the population
trial = lambda :  2500 <= sorted(randrange(10000) for i in range(5))[2] < 7500
print(sum(trial() for i in range(10000)) / 10000.0)    # compare with the beta distribution

# An exact permutation text
# Hypothesis testing
drug = [54, 73, 53, 70, 73, 68, 52, 65, 65]
placebo = [54, 51, 58, 44, 55, 52, 42, 47, 58, 46]
obs_diff = mean(drug) - mean(placebo)

# Null hypothesis is that there is no real difference between the drug
# and the placebo.  Which means that any observed difference was just
# due to chance and noise.

comb = drug + placebo

cases = 0
for i in range(10000):
    shuffle(comb)
    new_diff = mean(comb[:len(drug)]) - mean(comb[len(drug):])
    cases += new_diff >= obs_diff
print(cases / 10000, '<-- The p-value')
print('Accordlingly, we reject the null hypothesis')
print('and conclude the observed difference was not due to chance')