Simpleton

Lets proove Bayesian vs Frequentism

Sa 16 Juli 2016 — under ,

Bayesian vs Frequentist

There is a lot of ink spent on the topic by most prominent scholars of all time so I really dont want to add anything to the debate.

The 1 paper which influenced me in getting a practical overview of the topic and showing tangiblly was by Jake Vanderplass , where he explains with 2 examples

  • The table example
  • The coin example

It establishes from a practical prespective that Bayesian answer to the problem is near to the real world. But it made me thinking how do we know the real world answer? Either we need to play the table example for 1000 times or toss the coin 1000 times. The answer to the question is YES again which was provided by Jake himself in another (talk)[]. So in this work I try to experimentally validate both the experiments by running the examples 1000 times and finding the real values.

The Table example

Where frequentist and Bayesian differ. [http://www.rpgroup.caltech.edu/courses/PBoC_CSHL_2015/files_2015/articles/eddy%20bayesian%20methods.pdf]

Problem is to find the conditional probability of Bob wining given that Alice has won 5 and Bob 3. The game is won whoever reaches 6 points first!

P(Bob | A=5, B=3) ?

In [6]:
import scipy.stats as stats
import numpy as np
In [43]:
prob_bob_wining = np.array([])

for j in range(100):
    count_alice_win = 0
    count_bob_win = 0
    count_cond_hit = 0

    for i in range(10000):
        # Select random p 
        p = np.random.uniform(0,1)
        alice = 0
        bob = 0

        flag_condition = False
        
        # Loop till anyone wins. The game ends when anyone reaches 6 points
        while (alice < 6 and bob < 6):
            #The outcome of a single trial game
            outcome = np.random.binomial(1, p)
            
            #Alice wins if outcome is 1
            alice += outcome
            #Bob winsif outcome is 0
            bob += (1 - outcome)

            #The condition which is provided in the problem
            if (alice == 5 and bob == 3):
                flag_condition = True
                count_cond_hit += 1

        #Someone has won the game
        #checking if the condition was hit
        #in such condition how many times did bob win
        if flag_condition:
            if 6 == alice:
                count_alice_win += 1
            else:
                count_bob_win += 1
    
    #The divison of condition_hit to the time bob win is the
    #condition probability of Bob wining given P(Bob| A=5, B=3)
    prob_bob_wining = np.append(prob_bob_wining, count_bob_win/ count_cond_hit)

print ("P : ",p ," condition hit : ", count_cond_hit)
print ("Bob win after condition : ", count_bob_win)
print ("Alice win after condition : ", count_alice_win)
print ("Bob chance of wining ::   ", count_bob_win/count_cond_hit)
P :  0.15598397788391893  condition hit :  1170
Bob win after condition :  111
Alice win after condition :  1059
Bob chance of wining ::    0.09487179487179487
In [44]:
%matplotlib inline
import matplotlib.pyplot as plt
#Plotting the conditional probability of p(Bob| A=5, B=3) 
plt.hist(prob_bob_wining)
Out[44]:
(array([  1.,   4.,  11.,  13.,  20.,  20.,  21.,   7.,   2.,   1.]),
 array([ 0.07123034,  0.0753562 ,  0.07948206,  0.08360792,  0.08773378,
         0.09185964,  0.09598549,  0.10011135,  0.10423721,  0.10836307,
         0.11248893]),
 )

Bayesian view of Bob wining : Analytic Solution

{See the document for the proof}

This is a special case of Bayesian analysis in which we can can have a analytic solution. (Not true for all other cases of bayesian analysis)

$$ P( Bob | A=5, B=3 ) = \frac {B(5+a,3+b+3) }{B(5+a,3+b)} $$
In [3]:
import math
 
a_1 = 1 #prior
a_2 = 1 #prior for uniform beta distribution. Non informative prior

def beta_1(a,b):
     
    '''uses gamma function or inbuilt math.gamma() to compute values of beta function'''
     
    beta = math.gamma(a)*math.gamma(b)/math.gamma(a+b)
    return beta


print (beta_1(5 + a_1, 3+a_2+3) / beta_1(5 + a_1, 3 + a_2 ))
0.09090909090909091

Conclusion

Frequentist of Bob wining = $(5/3)^3 $ = 1/18 = 0.05

Bayesian of Bob wining = 0.09 $\sim$ 1/10

By Bootstrapping = 0.09-0.1 $\sim$ 1/10

In [ ]:
 

Coin example : Two heads

http://www.behind-the-enemy-lines.com/2008/01/are-you-bayesian-or-frequentist-or.html

Need to find

P(HH | H=10, T=4 ) ?

Bootstraping method

In [7]:
count_2_heads = 0
count_condition_hit = 0

percentage_of_2_heads = np.array([])
#for p in np.arange(0.1, 1, 0.1):
for j in range(100000):
    

    #for i in range(10000):
    #randomly find a p value
    p = np.random.uniform(0,1)

    #Toss the coint 14 times
    outcome = np.random.binomial(1,p, 14)

    #Check if the coin toss had 10 heads
    if (10 == outcome.sum() ):
        count_condition_hit += 1
        
        #Find the probability of 2 heads
        heads_2 = np.random.binomial(1,p, 2)
 
        # Check if there are 2 heads
        if (2 == heads_2.sum()):
            count_2_heads += 1

    #print ("For P : ", p ,"Condition Hit : ", count_condition_hit, "Count 2 heads : ", count_2_heads)
    #if (count_condition_hit):
    #    print ("Percentage of 2 Heads : ", (count_2_heads/ count_condition_hit))
    #else:
    #    print ("Precentage of 2 Heads : ", 0)
    if (count_condition_hit):
        percentage_of_2_heads = np.append(percentage_of_2_heads, (count_2_heads/ count_condition_hit))
    
print ("Percentage of 2 Heads : ", (count_2_heads/ count_condition_hit))
Percentage of 2 Heads :  0.47881743370923496

Bayesian Approach

$$ P ( HH | data ) = \frac {B(10+a+2,4+b) } {B(10+a,4+b)} $$

In [91]:
import math
 
def beta_1(a,b):
     
    '''uses gamma function or inbuilt math.gamma() to compute values of beta function'''
     
    beta = math.gamma(a)*math.gamma(b)/math.gamma(a+b)
    return beta


print (beta_1(10 + 1 + 2, 4+1) / beta_1(10 + 1 , 4 + 1 ))
0.4852941176470588

Conclusion

Frequentist answer of 2 Heads = $(10/14)^2 \sim 0.51 $

Bayesian answer of 2 Heads = 0.485

By Bootstrapping = 0.47


All content copyright 2014-2016 Deebul Nair unless otherwise noted. Licensed under Creative Commons.

Find me on Twitter, GitHub, email.