Wald Interval Lecture3 ICFP 2019.py
From Werner KRAUTH
Here is a simple Python2 program that checks the validity of the Wald interval for p=0.05 (corresponding to 1.96 sigma of the Gaussian distribution)
import math, random, pylab # # n: number of terms in the sum of Bernoulli variables # N_test: number of realizations of the sum of Bernoulli trials # theta: parameter of Bernoulli distribution # n_bin: Number of bins for theta # Here we use the value 1.96, appropriate for p=0.05 # n = 100 n_bins = 1000 N_test = 10000 x_values = [] y_values = [] for tt in range(2, n_bins - 1): theta = tt / float(n_bins) x_values.append(theta) Cover = 0.0 for iter in range(N_test): theta_hat = 0.0 for i in range(n): Upsilon = random.random() if Upsilon < theta: theta_hat += 1.0 / n if abs(theta_hat - theta) < 1.96 * math.sqrt(theta_hat * max(1 - theta_hat, 0.0) / n): Cover += 1.0 / N_test y_values.append(Cover) print theta, Cover, 'theta, Cover' pylab.title('Binomial proportion "standard 0.05 interval" $n = 100$') pylab.plot(x_values, y_values) pylab.xlabel('Bernoulli probability $ \\theta $') pylab.ylabel('Coverage probability $95$ \% standard interval') pylab.ylim(0.8,1.0) pylab.savefig('StandardInterval_test.png') pylab.show()