Comparación de las desigualdades de Bernstain y Chebyshev aplicadas a la distribución Bernoulli - la simulación en R da resultados inesperados

Question

Estoy tratando de comparar las desigualdades de Bernstein y Chebyshev aplicadas a la distribución Bernoulli con parámetro $p$ . Más específicamente - qué tan buenos son los límites que dan para diferentes tamaños de muestra. He escrito la simulación en R - los detalles se presentan a continuación. Los resultados y la pregunta precisa sobre ellos están en la parte inferior del post. Las desigualdades se establecen como sigue:

Dejemos que $\xi: Z \rightarrow \mathbb{R}$ sea una variable aleatoria con media $\mathbf{E}(\xi) = \mu$ y la varianza $\sigma^2(\xi) = \sigma^2 $ . Entonces, para cada $\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)}.$$ Si también $|\xi(z) - \mathbf{E}(\xi)| \leq M $ para casi todos los $z \in Z$ , entonces para cada $\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)} .$$

Escribí funciones que calculan el límite inferior de la probabilidad, que la diferencia entre la media empírica y la media no es mayor que $\varepsilon$ dado por ambas desigualdades:

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

El siguiente paso fue escribir la función que realiza un número de ensayos, donde cada ensayo realiza pasos:

1. generar m números extraídos según la distribución Bernoulli

2. calcular la diferencia entre la media empírica y la media

3. comprobar si la diferencia no supera $\varepsilon$

y luego calcula la fracción de ensayos en los que la diferencia entre la media empírica y la media no supera $\varepsilon$ . La función repite estos pasos 3 veces y también calcula los límites Chebyshev y Bernstein correspondientes.

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]))
}

Al final escribí una tabla, que recoge esas informaciones para diferentes tamaños de muestra. Aquí "muestra" significa vector poblado con tamaños de muestra que quiero comprobar.

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

Hago esta simulación varias veces con parámetros:

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

Cada vez obtuve resultados extraños. Este es uno de los resultados:

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

PROBLEMA: Las probabilidades empíricas son a veces menores que el límite inferior de probabilidad de Bernstain. En el ejemplo anterior ocurre para tamaños de muestra $m \geq 700$ . No sé por qué ocurre. He comprobado mis funciones y no encuentro ningún error. ¿Puede alguien ayudarme a resolver este problema o explicar la causa probable? ¿Quizás hay algún error numérico? ¿O el error es producido por la función 'rbinom'? He realizado una simulación similar para una distribución uniforme en $[0,1]$ y el resultado fue similar.

Answer 1

He ejecutado el código y no he visto ningún caso en el que se hayan superado los límites teóricos. Pero aun así, podría ocurrir, ya que el límite está en la probabilidad de que la media empírica se desvíe de la media real. Realiza el experimento $N$ veces, y debido a la mala suerte muchos de los medios empíricos podrían desviarse de $\mu$ De hecho, la proporción de juicios podría ser más de

$$2e^{ - \frac{m\varepsilon^2}{2(\sigma^2 + \frac{1}{3}M\varepsilon)}} $$

"violando" la desigualdad de Bernstein. Para un ejemplo sencillo, $P(\text{flipping 5 heads in a row}) \leq \frac{1}{10}$ (de hecho es $\frac{1}{32}$ ). Pero la realización de 1.000 ensayos de 5 volteos podría dar lugar a cualquier proporción de ensayos de 5 cabezas.