In [1]:
from __future__ import division

In [2]:
import random
def random_kid():
    return random.choice(['boy','girl'])

In [3]:
both_girls = 0
older_girl = 0
either_girl = 0

In [4]:
random.seed(0)

for _ in xrange(10000):
    younger = random_kid()
    older = random_kid()
    
    if older == 'girl':
        older_girl += 1
    if older == 'girl' and younger == 'girl':
        both_girls += 1
    if older == 'girl' or younger == 'girl':
        either_girl += 1

In [5]:
print 'P(both | older):', both_girls / older_girl
print 'P(both | either):', both_girls / either_girl


P(both | older): 0.514228456914
P(both | either): 0.341541328364

In [6]:
def uniform_pdf(x):
    return 1 if x >= 0 and x < 1 else 0

def uniform_cdf(x):
    if x < 0:
        return 0
    elif x < 1:
        return x
    else:
        return 1

Normal Distribution


In [8]:
import math
def normal_pdf(x, mu = 0, sigma = 1):
    sqrt_two_pi = math.sqrt(2 * math.pi)
    return (math.exp(-(x-mu)**2/ 2 / sigma**2) / (sqrt_two_pi * sigma))

def normal_cdf(x, mu = 0, sigma = 1):
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2

# binary search to make inverse normal
def inverse_normal_cdf(p, mu = 0, sigma = 1, tolerance = 0.0001):
    # make sure it's standard normal
    if mu != 0 or sigma != 1:
        return mu + sigma * inverse_normal_cdf(p, tolerance = tolerance)
    
    low_z, low_p = -10.0, 0
    hi_z, hi_p = 10.0, 1
    
    while hi_z - low_z > tolerance:
        mid_z = (low_z + hi_z) / 2
        mid_p = normal_cdf(mid_z)
        
        if mid_p < p:
            low_z, low_p = mid_z, mid_p
        elif mid_p > p:
            hi_z, hi_p = mid_z, mid_p
        else:
            break
    
    return mid_z

In [12]:
inverse_normal_cdf(0.05)


Out[12]:
-1.6448211669921875

Probability Distributions


In [13]:
def bernoulli_trial(p):
    return 1 if random.random() < p else 0

def binomial(n, p):
    return sum(bernoulli_trial(p) for _ in xrange(n))

In [16]:
# Investigate it
from collections import Counter

def make_hist(p, n, num_points):
    
    data = [binomial(n, p) for _ in xrange(num_points)]
    
    histogram = Counter(data)
    
    return histogram

In [17]:
make_hist(0.6, 5000, 100)


Out[17]:
Counter({2992: 4, 2990: 3, 3029: 3, 2957: 2, 2971: 2, 2973: 2, 2974: 2, 2978: 2, 2980: 2, 2982: 2, 2987: 2, 2997: 2, 2998: 2, 2999: 2, 3001: 2, 3013: 2, 3022: 2, 3027: 2, 3034: 2, 3069: 2, 2946: 1, 2947: 1, 3056: 1, 2953: 1, 2954: 1, 2956: 1, 2959: 1, 2960: 1, 2961: 1, 2962: 1, 2964: 1, 2965: 1, 2966: 1, 2967: 1, 2968: 1, 3097: 1, 3098: 1, 2975: 1, 2976: 1, 2977: 1, 2985: 1, 2986: 1, 2989: 1, 2994: 1, 3000: 1, 3002: 1, 3004: 1, 3005: 1, 3006: 1, 3061: 1, 3016: 1, 3017: 1, 3018: 1, 3023: 1, 3070: 1, 2900: 1, 3031: 1, 3032: 1, 3108: 1, 3037: 1, 3038: 1, 3039: 1, 3040: 1, 3041: 1, 3042: 1, 3044: 1, 3047: 1, 3052: 1, 3054: 1, 2928: 1, 2933: 1, 3062: 1, 3064: 1, 3065: 1, 2939: 1, 2942: 1})

In [ ]: