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