Comparison of Bernstain and Chebyshev inequalities applied to Bernoulli distribution - simulation in R gives unexpected results

I'm trying to compare Bernstein and Chebyshev inequalities applied to Bernoulli distribution with parameter $p$. More specifically - how good are bounds they give for different sample sizes. I wrote simulation in R - the details are presented below. The results and precise question about them are at the bottom of the post. The inequalities are stated as follows:

Let $\xi: Z \rightarrow \mathbb{R}$ be a random variable with mean $\mathbf{E}(\xi) = \mu$ and variance $\sigma^2(\xi) = \sigma^2$. Then, for every $\varepsilon>0$: $$\underset{\mathbf{z} \in Z^m}{\mathrm{Prob}} \Bigg\{ \Bigg| \frac{1}{m} \sum_{i=1}^{m} \xi(z_i) - \mu \Bigg| \leq \varepsilon \Bigg\} \geq 1 - \frac{\sigma^2}{m \varepsilon^2} \qquad \qquad \mbox{(Chebyshev)}.$$ If also $|\xi(z) - \mathbf{E}(\xi)| \leq M$ for almost all $z \in Z$, then for every $\varepsilon>0$ $$\underset{\mathbf{z} \in Z^m}{\mathrm{Prob}} \Bigg\{ \Bigg| \frac{1}{m} \sum_{i=1}^{m} \xi(z_i) - \mu \Bigg| \leq \varepsilon \Bigg\} \geq 1 - 2e^{ - \frac{m\varepsilon^2}{2(\sigma^2 + \frac{1}{3}M\varepsilon)}} \qquad \qquad \mbox{(Bernstein)} .$$

I wrote functions that calculate the lower bound on probability, that difference between empirical mean and mean is not greater than $\varepsilon$, given by both inequalities:

ChebyshevInequality <- function(m, epsilon, variance){
return( 1 - variance/(m*(epsilon^2)) )
}
BernsteinInequality <- function (m, epsilon, variance, M){
return( 1 - 2*exp( (-m*(epsilon^2))/(2*variance + 2*M*epsilon/3) ) )
}


Next step was to write function that performs a number of trials, where each trial preforms steps:

1. generate m numbers drawn accordingly to Bernoulli distribution

2. calculate the difference between empirical mean and mean

3. check if difference does not exceed $\varepsilon$

and then it calulates the fraction of trials, where difference between empirical mean and mean didn't exceed $\varepsilon$. The function repeats these steps 3 times and it also calculates corresponding Chebyshev and Bernstein bounds.

empiricalBernoulli <- function(trials, sample, p = .5, epsilon = .05)
{
m <- sample
mean <- p
var <- p*(1 - p)
M <- max(p, 1 - p)

does.not.exceed.epsilon <- c(logical(trials))
empirical <- c(numeric(3))

for(j in 1:3){

for(i in 1:trials){

observations <- rbinom(m, 1, p)
difference <- abs(sum(observations)/m - mean)
does.not.exceed.epsilon[i] <- (difference <= epsilon)

}

empirical[j] <- sum(does.not.exceed.epsilon)/trials
}

C <- ChebyshevInequality(m, epsilon, var)
B <- BernsteinInequality(m, epsilon, var, M)

return(data.frame(SampleSize = m, ChebyshevLowerBound = C, BernsteinLowerBound = B, Empirical1 = empirical[1], Empirical2 = empirical[2], empirical3 = empirical[3]))
}


At the end I wrote a table, that gathers those informations for different sample sizes. Here 'sample' means vector populated with sample sizes I want to check.

comparisonForBernoulli <- function(sample, trials = 1000, p = .5, epsilon = .05){

k <- length(sample)

df <- data.frame(SampleSize = numeric(k), ChebyshevLowerBound = numeric(k), BernsteinLowerBound = numeric(k), Empirical1 = numeric(k), Empirical2  = numeric(k), Empirical3  = numeric(k))

for(i in 1:k){
df[i,] <- empiricalBernoulli(trials, sample[i], p , epsilon)
}
return(df)
}


I run this simulation several times with parameters:

sample1 <- seq(100,1000,100)
comparisonForBernoulli(sample = sample1, trials = 1000, p = .5, epsilon = .05)


Each time I obtained weird results. This is one of the outputs:

SampleSize ChebyshevLowerBound BernsteinLowerBound Empirical1 Empirical2 Empirical3
1         100        2.220446e-16          -0.2327855      0.709      0.685      0.679
2         200        5.000000e-01           0.2401200      0.850      0.823      0.854
3         300        6.666667e-01           0.5316155      0.910      0.922      0.913
4         400        7.500000e-01           0.7112912      0.942      0.950      0.957
5         500        8.000000e-01           0.8220420      0.974      0.973      0.967
6         600        8.333333e-01           0.8903080      0.991      0.982      0.989
7         700        8.571429e-01           0.9323866      0.993      0.990      0.994
8         800        8.750000e-01           0.9583236      0.995      0.993      0.999
9         900        8.888889e-01           0.9743110      0.997      0.999      0.995
10       1000        9.000000e-01           0.9841655      0.999      0.997      0.999


PROBLEM: The empirical probabilities are sometimes smaller than Bernstain's lower bound on probability. In the example above it happens for sample sizes $m \geq 700$. I do not know why does it happen. I checked my functions and I can't find any mistakes. Can someone help me to resolve this issue or explain the probable cause? Maybe there are some numeric errors? Or the error is produced by function 'rbinom'? I run similar simulation for uniform distibution on $[0,1]$ and the output was similar.

• I haven't the time to trace through your code, but a glance at it suggests you might be confusing the true mean and SD of the underlying distributions with their estimates from the simulated samples.
– whuber
Commented Jul 1, 2015 at 13:46
• Totally unrelated R comment: c(numeric(100)) etc. is repeating similar thing twice, because it says: "make a vector of a numeric vector of length 100" so using numeric(100) (more efficient if you know the length in advance) or simply c() (empty vector of length 0) or x <- NULL (if you do not want to assume anything about the vector) is enough.
– Tim
Commented Jul 2, 2015 at 4:40

I ran the code and didn't see any cases for which the theoretical bounds were breached. But even so, it could happen, as the bound is on the probability that the empirical mean deviates from the true mean. Run the experiment $N$ times, and due to bad luck many of the empirical means could deviate from $\mu$, in fact, the proportion of trials could be more than
$$2e^{ - \frac{m\varepsilon^2}{2(\sigma^2 + \frac{1}{3}M\varepsilon)}}$$
"violating" Bernstein's inequality. For a simple example, $P(\text{flipping 5 heads in a row}) \leq \frac{1}{10}$ (it is in fact $\frac{1}{32}$). But performing 1000 trials of 5 flips could result in any proportion of 5-head trials.