Last active
December 13, 2023 18:20
-
-
Save endolith/9e73ea30511f0befb9b336526f34d0a6 to your computer and use it in GitHub Desktop.
Confirm the formula for estimating the variance of an estimated probability. (Trying to understand error bars for https://github.com/endolith/elsim/issues/13)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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 | |
""" |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment