In [1]:

    

alpha <- 0.05
p <- 0.4
sampleSize <- 100000


    

In [2]:

    

startExp <- 0
endExp <- 4


    

In [3]:

    

IsInInterval <- function(center, width, x) {
    # Return TRUE if x is in interval (center-width, center+width)
    return ((center - width) < x && x < (center + width))
}


    

In [4]:

    

Epsilon <- function(alpha, n) {
    return (sqrt(1/(2*n) * log(2/alpha)))
}


    

In [5]:

    

CntBernoulliOverlap <- function(n, p, sampleSize, alpha) {
    # Return how often the Bernoulli-parameter p is in interval
    #
    # Args:
    #   n: Bernoulli Trials
    #   p: Bernoulli-parameter
    #   sampleSize:
    #   alpha:
    eps <- Epsilon(alpha, n)
    
    sampleBernSum <- rbinom(sampleSize, n, p)
    sampleMean <- sampleBernSum / n

    sampleInInterval <- sapply(sampleMean, IsInInterval, width=eps, x=p)
    return (sum(sampleInInterval))
}


    

In [6]:

    

ns <- 10^(startExp:endExp)
ns


    

    
Out[6]:





	
1
	
10
	
100
	
1000
	
10000

In [7]:

    

counts <- sapply(ns, CntBernoulliOverlap, p=p, sampleSize=sampleSize, alpha=alpha) 
counts


    

    
Out[7]:





	
100000
	
99845
	
99446
	
99395
	
99454

In [8]:

    

plot(ns, counts, log="x", type="h")


    

In [9]:

    

epsilons <- sapply(ns, Epsilon, alpha=alpha)
lengths = 2 * epsilons
lengths


    

    
Out[9]:





	
2.71620303148124
	
0.858938816693475
	
0.271620303148124
	
0.0858938816693475
	
0.0271620303148124

In [10]:

    

plot(ns, lengths, log="x")


    

In [ ]: