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:
generate m numbers drawn accordingly to Bernoulli distribution
calculate the difference between empirical mean and mean
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.
c(numeric(100))
etc. is repeating similar thing twice, because it says: "make a vector of a numeric vector of length 100" so usingnumeric(100)
(more efficient if you know the length in advance) or simplyc()
(empty vector of length 0) orx <- NULL
(if you do not want to assume anything about the vector) is enough. $\endgroup$