Confirm the formula for estimating the variance of an estimated probability. (Trying to understand error bars for https://github.com/endolith/elsim/issues/13)

estimated_variance_of_estimated_probability.py

"""
Confirm the formula for estimating the variance of an estimated probability.

Created on Fri May 27 11:06:43 2022
"""
import numpy as np


def func(x):
    """
    Return True with an unknown probability
    """
    if x < 0.35098271064435571344:  # From random.org lol
        return True
    else:
        return False


def var(p, n):
    """
    Estimate the variance of an estimate.

    Parameters
    ----------
    p : float
        Proportion of trials that resulted in success.
    n : int
        The number of trials performed.

    Returns
    -------
    v : float
        The estimated variance of the estimate.

    """
    r = 1.0  # TODO: Make this variable for all
    return p*(r - p)/n


def monte_carlo(func, n):
    """
    Run `n` trials of `func` and estimate the probability it returns True.

    Parameters
    ----------
    func : function
        A function that takes an input in the range [0.0, 1.0) and returns True
        or False.
    n : int
        The number of trials to perform.

    Returns
    -------
    p : float
        The proportion of trials that returned True.

    """
    hits = 0
    inputs = np.random.random_sample(n)
    for x in inputs:
        if func(x):
            hits += 1

    return hits / n

n = 100_000
m = 10_000
estimates = []
for i in range(m):
    p = monte_carlo(func, n)
    estimates.append(p)
    if i < 10:
        print(f'Estimated probability {p:.5f} with estimated variance '
              f'{var(p, n):.3e}')
    elif i == 10:
        print('[More trials]...')

print(f'\nActual variance of outputs of all trials: {np.var(estimates):.3e}')

"""
Output with n = 100_000, m = 10_000:

Estimated probability 0.35028 with estimated variance 2.276e-06
Estimated probability 0.35364 with estimated variance 2.286e-06
Estimated probability 0.34948 with estimated variance 2.273e-06
Estimated probability 0.35267 with estimated variance 2.283e-06
Estimated probability 0.34873 with estimated variance 2.271e-06
Estimated probability 0.34995 with estimated variance 2.275e-06
Estimated probability 0.35047 with estimated variance 2.276e-06
Estimated probability 0.35072 with estimated variance 2.277e-06
Estimated probability 0.35314 with estimated variance 2.284e-06
Estimated probability 0.35018 with estimated variance 2.276e-06
[More trials]...

Actual variance of outputs of all trials: 2.280e-06
"""