# Wald Interval Lecture3 ICFP 2019.py

Revision as of 09:42, 24 September 2019; view current revision

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()