En este ejemplo veremos cómo podemos calcular los intervalos de confianza vistos en clase a partir de una muestra. Para ello usaremos el conjunto de datos mtcars de R y no centraremos en la variable mpg que representa las millas por galon correspondientes a cada una de las 32 observaciones.

Supongamos que se sabe que la variable sigue una distribución Normal y queremos calcular un intervalo del 95% de confianza para la media. Como no tenemos información sobre la varianza, el intervalo que saremos es el siguiente: \[ \bar{x} \pm \frac{S}{\sqrt{n}}t_{n-1}^{1-\alpha/2} \] Primero cargamos los datos y calculamos una a una las cantidades que necesitamos:

data("mtcars")
media<- mean(mtcars$mpg)
print(media)
## [1] 20.09062
n <- length(mtcars$mpg)
error.estandar<- sd(mtcars$mpg)/sqrt(n)
alpha = 0.05
grados.lib <- n - 1
cuantil.t = qt(1-(alpha/2), grados.lib)

Recordemos que a \(\frac{S}{\sqrt{n}}t_{n-1}^{1-\alpha/2}\) se le conoce como margen de error, primero calculamos el margen de error y finalmente el intervalo:

margen.error <- cuantil.t * error.estandar
print(margen.error)
## [1] 2.172946
limite.inferior <- media - margen.error
limite.superior<- media + margen.error
print(c(limite.inferior,limite.superior))
## [1] 17.91768 22.26357

A continuación haremos un ejemplo simulado para observar la diferencia que provoca conocer o no el valor de \(\sigma\).

Primero simulamos una muestra de tamaño 30 proveniente de una población con distribución normal de media 5 y varianza 9

n<-30
mu<-5
sigma<-3
muestra<-rnorm(n,mu,sigma)

Si suponemos que \(\sigma=3\) es conocida, entonces para encontrar un intervalo del 95% de confianza para la media tendríamos que usar \[ \bar{x} \pm z_{0.975}\frac{\sigma}{\sqrt{n}} \] ya que \(\alpha=0.05\). Usando la muestra obtenemos lo siguiente:

margen.error<-qnorm(0.975,0,1)*(sigma/sqrt(n))
intervalo1<-c(mean(muestra)-margen.error,mean(muestra)+margen.error)
intervalo1
## [1] 3.650909 5.797942

Ahora, si suponemos que \(\sigma\) es desconocida, el intervalo sería

margen.error2<-qt(0.975,15)*(sd(muestra)/sqrt(n))
intervalo2<-c(mean(muestra)-margen.error2,mean(muestra)+margen.error2)
intervalo2
## [1] 3.730427 5.718424

Haga diferentes pruebas cambiando los valores de \(n\), \(\mu\) y \(\sigma\) para observar qué pasa con los intervalos. ¿Hay diferencia en los intervalos cuando se conoce o no el valor de \(\sigma\)?

A continuación calcularemos un intervalo del 95% de confianza para \(\sigma\). En este caso, es común que se utilice un intervalo simétrico, es decir, se divide al nivel de significancia \(\alpha\) en dos partes iguales (como si la distribución de la cantidad pivotal fuera simétrico) y se obtiene que el intervalo es el siguiente: \[ \left( \sqrt{\frac{(n-1)S^2}{\chi^2_{0.975, n-1}}}, \sqrt{\frac{(n-1)S^2}{\chi^2_{0.025, n-1}}} \right). \]

Usando la muestra se obtiene el valor del intervalo:

qc=qchisq(c(.975,.025),(n-1))
intervalo3<-sqrt(((n-1)*(sd(muestra)^2))/qc)
intervalo3
## [1] 2.03426 3.43378

Sin embargo, recordemos que los valores \(\chi^2_{0.975, n-1}\) y \(\chi^2_{0.025, n-1}\) son los cuantiles \(0.025\) y \(0.975\) que se eligen para acumular una probabilidad de \(1-\alpha=0.95\). Entonces tambien se pudieron haber utilizado los cuantiles \(\chi^2_{0.98, n-1}\) y \(\chi^2_{0.03, n-1}\) para acumular el \(0.95\), en tal caso se hubiera obtenido el siguiente intervalo:

qc=qchisq(c(.98,.03),(n-1))
intervalo4<-sqrt(((n-1)*(sd(muestra)^2))/qc)
intervalo4
## [1] 2.013010 3.391078

La diferencia principal entre los intervalos anteriores es la longitud:

intervalo3[2]- intervalo3[1]
## [1] 1.39952
intervalo4[2]- intervalo4[1]
## [1] 1.378068

¿Qué podemoha hacer para obtener el intervalo de longitud mínima?. Sabemos que cuando la distribución de Q es simétrica, el intervalo de longitud mínima se obtiene al divider \(\alpha\) en dos partes iguales, sim embargo en este caso es necesario encontrar la forma óptima de dividir \(\alpha\). Hay mas de una forma de encontrar que valores usar, aquí tomaremos varios valores entre 0 y 0.05, calcularemos los intervalos y veremos cual de ellos tiene longitud mínima:

lp = seq(0.001, .049, by=.001) #dividimos el intervalo (0,0.5) en partes iguales
m = length(lp) #Número de intervalos que calcularemos
len=numeric(m) #Aquí iremos guardando las longitudes
for(i in 1:m) {
  L = qchisq(lp[i], (n-1)) #Primer cuantil
  U = qchisq(.95+lp[i], (n-1)) #segundo cuantil
  lcl = sqrt((n-1)*sd(muestra)^2/U) #limite inferior del intervalo
  ucl = sqrt((n-1)*sd(muestra)^2/L)  #limite superior del intervalo
  len[i] = ucl-lcl } #longitud
plot(lp, len, type="l", lwd=2) #graficamos las longitudes en función del valor usado

min(len) #Longitud mínima
## [1] 1.366464
alpha<-lp[len==min(len)] #Vemos qué valor nos dio la longitud mínima
alpha
## [1] 0.037
qc=qchisq(c(alpha+0.95,alpha),(n-1))
sqrt(((n-1)*(sd(muestra)^2))/qc)#el intervalo
## [1] 1.974914 3.341378