difference_of_biased_estimators.py

# Michael A. Alcorn (airalcorn2@gmail.com)

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import random

fixed_rates = {"p_1": 0.3, "p_2": 0.5}

find_rates = {"biased": {"fixed": 0.7, "non_fixed": 0.9},
              "unbiased": {"fixed": 0.9, "non_fixed": 0.9}}

trials = 10000
sampling_attempts = 1000
data = []
for trial in range(trials):
    
    counts = {}
    for bias in find_rates:
        counts[bias] = {}
        for p in fixed_rates:
            counts[bias][p] = {"fixed": 0, "non_fixed": 0}
    
    for bias in find_rates:
        for sample in range(sampling_attempts):
            for p in fixed_rates:
                if random.random() < fixed_rates[p]:
                    # Fixed locus.
                    if random.random() < find_rates[bias]["fixed"]:
                        counts[bias][p]["fixed"] += 1
                else:
                    # Non-fixed locus.
                    if random.random() < find_rates[bias]["non_fixed"]:
                        counts[bias][p]["non_fixed"] += 1
    
    row = {}
    for bias in find_rates:
        
        (p_1_fixed, p_1_non_fixed) = (counts[bias]["p_1"]["fixed"], counts[bias]["p_1"]["non_fixed"])
        p_1_prop = p_1_fixed / (p_1_fixed + p_1_non_fixed)
        
        (p_2_fixed, p_2_non_fixed) = (counts[bias]["p_2"]["fixed"], counts[bias]["p_2"]["non_fixed"])
        p_2_prop = p_2_fixed / (p_2_fixed + p_2_non_fixed)
        
        row[bias] = p_1_prop - p_2_prop
    
    data.append(row)

df = pd.DataFrame(data)
sns.distplot(df["biased"], label = "biased")
sns.distplot(df["unbiased"], label = "unbiased")
plt.legend()
plt.show()