Temas

Esta sección esta basada en el libro An Introduction to the Bootstrap de Effron y Tibshirani

Muestras aleatorias

Supongamos que tenemos una población finita o universo \(U\), conformado por unidades individuales \(U_1,U_2,...,U_N\) cada una tiene la misma probabilidad de ser seleccionada en una extracción aleatoria. Las unidades individuales \(U_i\) tienen propiedades que nos gustaría aprender (opinión política,…). Debido a que es muy difícil y caro examinar cada unidad en \(U\) seleccionamos una muestra aleatoria.

Una muestra aleatoria de tamaño \(n\) se define como una colección de \(n\) unidades \(u_1,...,u_n\) seleccionadas aleatoriamente de un universo \(U\).

En principio el proceso de muestreo es como sigue:

  1. Seleccionamos \(n\) enteros de manera independiente \(j_1,...,j_n\) (con probabilidad \(1/N\)), cada uno de ellos asociado a un número entre \(1\) y \(N\).

  2. Los enteros determinan las unidades que seleccionamos: \(u_1=U_{j_1},u_2=U_{j_2},...,u_n=U_{j_n}\).

En la práctica el proceso de selección suele ser más complicado y la población definición de la población \(U\) suele ser deficiente; sin embargo el marco conceptual sigue siendo útil para entender la inferencia estadística.

Observación: Nuestra definición de muestra aleatoria permite que una unidad particular \(U_i\) aparezca más de una vez, podríamos evita esto si realizamos un muestreo sin remplazo; sin embargo, es un poco más sencillo permitir repeticiones y si el tamaño de la muestra \(n\) es mucho más chico que la población \(N\), la probabilidad de muestrear la misma unidad más de una vez es chica.

Una vez que selecciona una muestra aleatoria \(u_1,...,u_n\) obtenemos una o más medidas de interés para cada unidad. Los datos observados son la colección de medidas \(x_1,...,x_n\), que también denotaremos \(\textbf{x} = (x_1,...,x_n)\).

Podemos imaginar también, obtener las medias de interés de cada unidad en la población \(U_1,U_2,...,U_N\), obteniendo así los valores \(X_1,...,X_N\), esto sería un censo, y denotamos al conjunto de mediciones de la población por \(\mathcal{X}\). El objetivo de la inferencia estadística es expresar lo que hemos aprendido de la población \(\mathcal{X}\) a partir de los datos observados \(\textbf{x}\). En particular, vamos a usar el bootstrap para determinar que tan preciso es un estadístico (e.g. media o mediana) calculado de la muestra \(x_1,...,x_n\) estima la cantidad correspondiente en la población.

Veamos un ejemplo artificial donde tenemos una muestra de 500 escuelas primarias del DF tomada de un universo de 3,311 escuelas,

set.seed(16021)
n <- 500
prim_muestra <- sample_n(prim, n, replace = TRUE)
head(prim_muestra)
## Source: local data frame [6 x 6]
## 
##        clave      turno       tipo mun   esp3   esp6
## 1 09DPR3043G   MATUTINO    GENERAL  16 538.89 542.61
## 2 09PPR1369C   MATUTINO PARTICULAR  18 670.16 591.52
## 3 09DBN0026X   NOCTURNO    GENERAL  26 629.82   0.00
## 4 09DPR2982T VESPERTINO    GENERAL  16 536.35 563.03
## 5 09DPR2906N   MATUTINO    GENERAL  26 539.10 588.37
## 6 09PPR1515X   MATUTINO PARTICULAR  18 542.62 617.60

para cada escuela en la muestra tenemos la medida \(x_i\), conformada por el promedio de las calificaciones en español de los alumnos de tercero y sexto de primaria (prueba ENLACE 2010): \[x_i=(esp3_i, esp6_i)\]

Este ejemplo es artificial pues contamos con un censo de las escuelas, sin embargo es común contar únicamente con la muestra, esta tiene una media de 571.6, con un error estándar estimado de 4.01. Debido a que nuestro ejemplo es artificial podemos comparar con la población, la media de las 3311 escuela es 574.8.

El principio del plug-in

En ocasiones, los problemas de inferencia estadística involucran la estimación de algún aspecto de una distribución de de probabilidad \(P\) en base a una muestra aleatoria obtenida de \(P\). La función de distribución empírica \(P_n\) es una estimación de la distribución completa \(P\), por lo que una manera inmediata de estimar aspectos de \(P\) (e.g media o mediana) es calcular el aspecto correspondiente de \(P_n\).

Podemos comparar el histograma de la distribución completa con el histograma de la distribución empírica para el ejemplo de las calificaciones de la prueba ENLACE.

claves_muestra <- prim_muestra$clave
prim_tidy <- prim %>% 
  gather(grado, calif, esp3, esp6) %>%
  mutate(muestra = ifelse(clave %in% claves_muestra, "muestra", "población"))

prim_plot <- prim_muestra  %>% 
  gather(grado, calif, esp3, esp6) %>%
  mutate(muestra = "población") %>%
  rbind(prim_tidy)

ggplot(prim_plot, aes(x = calif)) +
  geom_histogram(aes(y = ..density..), binwidth = 20, fill = "darkgray") +
  facet_grid(grado ~ muestra)

Cuando la variable de interés toma pocos valores es fácil ver la distribución empírica: supongamos que la medición de las unidades que nos interesa es la variable tipo de escuela, entonces la distribución empírica en la muestra es

table(prim_muestra$tipo) / n
## 
##    GENERAL PARTICULAR 
##      0.654      0.346

Vale la pena notar que pasar de la muestra desagregada a la distribución empírica (lista de valores y la proporción que ocurre cada una en la muestra) no conlleva ninguna pérdida de información: el vector de frecuencias observadas es un estadístico suficiente para la verdadera distribución. Esto quiere decir que toda la información de \(P\) contenida en el vector de observaciones \(\textbf{x}\) está también contenida en \(P_n\).

Nota: el teorema de suficiencia asume que las observaciones \(\textbf{x}\) son una muestra aleatoria de la distribución \(P\), este no es siempre el caso (e.g. si tenemos una serie de tiempo).

Cuando aplicamos teoría estadística a problemas reales, es común que las respuestas estén dadas en términos de distribuciones de probabilidad. Por ejemplo, podemos preguntarnos que tan correlacionados están los resultados de las pruebas de español correspondientes a 3^o y 6^o. Si conocemos la distribución de probabilidad \(P\) contestar esta pregunta es simplemente cuestión de aritmética, el coeficiente de correlación poblacional esta dado por:

\[corr(y,z) = \frac{\sum_{j=1}^{N}(Y_j - \mu_y)(Z_j-\mu_z)} {[\sum_{j=1}^{N}(Y_j - \mu_y)^2\sum_{j=1}^{N}(Z_j - \mu_z)^2]^{1/2}}\]

en nuestro ejemplo \((Y_j,Z_j)\) son el j-ésimo punto en la población de escuelas primarias \(\mathcal{X}\), \(\mu_y=\sum Y_j/3311\) y \(\mu_z=\sum Z_j/3311\).

ggplot(prim, aes(x = esp3, y = esp6)) +
  geom_point(alpha = 0.5)

cor(prim$esp3, prim$esp6)
## [1] 0.2396605

Si no tenemos un censo debemos inferir, podríamos estimar la correlación \(corr(y,z)\) a través del coeficiente de correlación muestral: \[\hat{corr}(y,z) = \frac{\sum_{j=1}^{n}(y_j - \hat{\mu}_y)(z_j-\hat{\mu}_z)} {[\sum_{j=1}^{n}(y_j - \hat{\mu}_y)^2\sum_{j=1}^{n}(z_j - \hat{\mu}_z)^2]^{1/2}}\]

cor(prim_muestra$esp3, prim_muestra$esp6)
## [1] 0.2798717

Otros ejemplos de estimaciones plug-in:

median(prim_muestra$esp3)
## [1] 561.855

\[\theta=\frac{1}{N}\sum_{j=1}^N I_{\{Y_i>700\}}\]

donde \(I_{\{\cdot\}}\) es la función indicadora.

Hacemos la estimación plug-in \(\hat{\theta}\):

sum(prim_muestra$esp3 > 700) / n
## [1] 0.046

Volvamos ahora al ejemplo de los 100 lanzamientos de un dado, en este ejemplo obteníamos la siguiente distribución empírica:

dado <- read.table("../04-Probabilidad/data/dado.csv", header=TRUE, quote="\"")
table(dado) / n
## dado
##     1     2     3     4     5     6 
## 0.026 0.038 0.020 0.034 0.028 0.054

En este caso no tenemos un censo, solo contamos con la muestra. Una pregunta de inferencia que surge de manera natural es si el dado es justo, esto es, si la distribución que generó esta muestra tiene una distribución \(P = (1/6, 1/6, 1/6,1/6, 1/6, 1/6)\). Para resolver esta pregunta, debemos hacer inferencia de la distribución empírica.

Antes de proseguir repasemos dos conceptos importantes: parámetros y estadísticos:

Un parámetro es una función de la distribución de probabilidad \(\theta=t(P)\), mientras que un estadístico es una función de la muestra \(\textbf{x}\).

Por ejemplo, la \(corr(y,z)\) es un parámetro de \(P\) y \(\hat{corr}(x,y)\) es un estadístico basado en \(\textbf{x}\).

Entonces:

El principio del plug-in es un método para estimar parámetros a partir de muestras; la estimación plug-in de un parámetro \(\theta=t(P)\) se define como: \[\hat{\theta}=t(P_n).\]

Una pregunta natural es: ¿qué tan bien funciona el principio del plug-in? Usualmente es muy bueno cuando la única información disponible de \(P\) es la muestra \(\textbf{x}\), bajo esta circunstancia \(\hat{\theta}=t(P_n)\) no puede ser superado como estimador de \(\theta=t(P)\), al menos no en el sentido asintótico de teoría estadística \((n\to\infty)\).

El principio del plug-in provee de una estimación más no habla de precisión, por ello usaremos el bootstrap para estudiar el sesgo y el error estándar del estimador plug-in \(\hat{\theta}=t(P_n)\), la maravilla del bootstrap es que produce errores estándar y sesgos de manera automática, sin importar que tan complicada es la función \(t(P)\).

Errores estándar y sus estimaciones

Los estadísticos como \(\hat{\theta}=t(P_n)\) suelen ser el primer paso en el análisis de datos, el siguiente paso es investigar la precisión de las estimaciones; el bootstrap es un método para calcular precisión de estimaciones que se vale del principio del plug-in para estimar el error estándar de una estadística.

Ejemplo: el error estándar de una media

Supongamos que \(x\) es una variable aleatoria que toma valores en los reales con distribución de probabilidad P. Denotamos por \(\mu_P\) y \(\sigma_P^2\) la media y varianza de P,

\[\mu_P = E_P(x),\] \[\sigma_P^2=var_P(x)=E_P[(x-\mu_P)^2]\]

en la notación enfatizamos la dependencia de la media y varianza en la distribución \(P\).

Ahora, sea \((x_1,...,x_n)\) una muestra aleatoria de \(P\), de tamaño \(n\), la media de la muestra \(\bar{x}=\sum_{i=1}^nx_i/n\) tiene esperanza \(\mu_P\) y varianza \(\sigma_P^2/n\).

En palabras: la esperanza de \(\bar{x}\) es la misma que la esperanza de \(x\), pero la varianza de \(\bar{x}\) es \(1/n\) veces la varianza de \(x\), así que entre mayor es la \(n\) tenemos una mejor estimación de \(\mu_P\).

El error estándar denota la desviación estándar de una estadística. En el caso de la media \(\bar{x}\), el error estándar, que denotamos \(se_P(\bar{x})\), es la raíz de la varianza de \(\bar{x}\), \[se_P(\bar{x}) = [var_P(\bar{x})]^{1/2}= \sigma_P/ \sqrt{n}.\]

En este punto podemos usar el principio del plug-in, simplemente sustituimos \(P_n\) por \(P\) y obtenemos, primero, una estimación de \(\sigma_P\): \[\hat{\sigma}=\hat{\sigma}_{P_n} = \bigg\{\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2\bigg\}^{1/2}\]

de donde se sigue la estimación del error estándar: \[\hat{se}(\bar{x})=\sigma_{P_n}/\sqrt{n}=\bigg\{\frac{1}{n^2}\sum_{i=1}^n(x_i-\bar{x})^2\bigg\}^{1/2}\]

Notemos que usamos el principio del plug-in en dos ocasiones, primero para estimar la esperanza \(\mu_P\) mediante \(\mu_{P_n}\) y luego para estimar el error estándar \(se_P(\bar{x})\) mediante \(se_{P_n}(\bar{x})\). En el caso de la media \(\hat{\theta}=\bar{x}\) la aplicación del principio del plug-in para el cálculo de errores estándar es inmediata; sin embargo, hay estadísticas para las cuáles no es fácil aplicar este método y es ahí cuando aplicaremos el bootstrap.

Antes de pasar al bootstrap podemos preguntarnos: ¿porqué tanto énfasis en el error estándar? El error estándar es la manera más común para describir la precisión de una estadística. En términos generales, esperamos que \(\bar{x}\) este a una distancia de \(\mu_P\) menor a un error estándar el 68% del tiempo, y a menos de 2 errores estándar el 95% del tiempo. Estos porcentajes están basados el teorema central del límite que nos dice que bajo ciertas condiciones (bastante generales) de \(P\) la distribución de \(\bar{x}\) se aproximará a una distribución normal: \[\bar{x} \overset{\cdot}{\sim} N(\mu_P,\sigma_P^2/n)\]

El estimador bootstrap del error estándar

Supongamos que tenemos una muestra aleatoria \(\textbf{x}=(x_1,x_2,...,x_n)\) proveniente de una distribución de probabilidad desconocida \(P_n\) y deseamos estimar un parámetro \(\theta = t(P)\) con base en la muestra. Para esto, calculamos una estimación \(\hat{\theta}=s(\textbf{x})\) (la estimación puede ser la estimación plug-in \(t(P_n)\) pero también puede ser otra). Entonces podemos usar bootstrap para calcular el error estándar de la estimación.

Definimos una muestra bootstrap como una muestra aleatoria de tamaño \(n\) que se obtiene de la distribución empírica \(P_n\) y la denotamos \[\textbf{x}^* = (x_1^*,...,x_n^*).\]

La notación de estrella indica que \(\textbf{x}^*\) no son los datos \(\textbf{x}\) sino una versión de remuestreo de \(\textbf{x}\).

Otra manera de frasearlo: Los datos bootsrtap \(x_1^*,...,x_n^*\) son una muestra aleatoria de tamaño \(n\) seleccionada con reemplazo de la población de \(n\) objetos \((x_1,...,x_n)\).

Ahora, a cada muestra bootstrap \(\textbf{x}^*\) le corresponde una replicación \(\hat{\theta}^*=s(\textbf{x}^*).\)

La estimación bootstrap de \(se_{P}(\hat{\theta}^*)\), esto es, el error estándar de un estadístico \(\hat{\theta}\) es una estimación plug-in en donde la distribución empírica \(P_n\) toma el lugar de la distribución desconocida \(P\): el estimador bootstrap de \(se_P(\hat{\theta})\) se define como: \[se_{P_n}(\hat{\theta}^*)\] en otras palabras, la estimación bootstrap de \(se_P(\hat{\theta})\) es el error estándar de \(\hat{\theta}\) para conjuntos de datos de tamaño \(n\) seleccionados de manera aleatoria de \(P_n\).

La fórmula \(se_{P_n}(\hat{\theta}^*)\) no existe para casi ninguna estimación que diferente de la media, por lo que recurrimos a la técnica computacional bootstrap: el algoritmo funciona seleccionando distintas muestras bootstrap, evaluando la replicación bootstrap correspondiente y estimando el error estándar de \(\hat{\theta}\) mediante la desviación estándar empírica de las replicaciones. El resultado es la estimación bootstrap del error estándar, que denotamos \(\hat{se}_B\), donde \(B\) es el número de muestras bootstrap usadas.

Algoritmo bootstrap para estimar errores estándar

  1. Selecciona \(B\) muestras bootsrtap independientes: \[\textbf{x}^{*1},..., \textbf{x}^{*B}\].
  2. Evalúa la replicación bootstrap correspondiente a cada muestra bootstrap: \[\hat{\theta}^{*b}=s(\textbf{x}^{*b})\] para \(b=1,2,...,B.\)
  3. Estima el error estándar \(se_P(\hat{\theta})\) usando la desviación estándar muestral de las \(B\) replicaciones: \[\hat{se}_B = \bigg\{\frac{\sum_{b=1}^B[\hat{\theta}^*(b)-\hat{\theta}^*(\cdot)]^2 }{B-1}\bigg\}^{1/2}\]
donde $^*()=_{b=1}^B }^*(b)/B $.

Conforme el número de replicaciones \(B\) aumenta \[\hat{se}_B\approx se_{P_n}(\hat{\theta})\] este hecho equivale a decir que la desviación estándar empírica se acerca a la desviación estándar poblacional conforme crece el número de muestras. La población en este caso es la población de valores \(\hat{\theta}^*=s(x^*)\)

Escribimos una función para calcular el error estándar de una media usando replicaciones bootstrap:

mediaBoot <- function(x){ 
  # x: variable de interés
  # n: número de replicaciones bootstrap
  n <- length(x)
  muestra_boot <- sample(x, size = n, replace = TRUE)
  mean(muestra_boot) # replicacion bootstrap de theta_gorro
}
thetas_boot <- rdply(1000, mediaBoot(prim_muestra$esp3))
sd(thetas_boot$V1)
## [1] 4.025791

y se compara con

se <- function(x) sqrt(sum((x - mean(x)) ^ 2)) / length(x)
se(prim_muestra$esp3)
## [1] 4.137324

Considera el coeficiente de correlación muestral entre la calificación de \(y=\)español 3 y la de \(z=\)español 6: \(\hat{corr}(y,z)=0.9\). ¿Qué tan preciso es esta estimación?

¿Cuántas replicaciones bootstrap (B)?

La estimación bootstrap ideal es un resultado asintótico \(B=\infty\), en esta caso \(\hat{se}_B\) iguala la estimación plug-in \(se_{P_n}\). En la práctica para elegir el tamaño de \(B\) debemos considerar que buscamos las mismas propiedades para la estimación de un error esrándar que para cualquier estimación: poco sesgo y desviación estándar chica. El sesgo de la estimación bootstrap del error estándar suele ser bajo y el error estándar está

Reglas de dedo (Effron y Tibshirani):

  1. Incluso un número chico de replicaciones bootstrap, digamos \(B=25\) es informativo, y \(B=50\) con frecuencia es suficiente para dar una buena estimación de \(se_P(\hat{\theta})\).

  2. En pocos casos es necesario realizar más de \(B=200\) replicaciones cuando se busca estimar error estándar.

seMediaBoot <- function(x, B){
  thetas_boot <- rdply(B, mediaBoot(x))
  sd(thetas_boot$V1)
}

B_muestras <- c(5, 25, 50, 100, 200, 400, 800, 1600, 3200)
sapply(B_muestras, function(i) seMediaBoot(x = prim_muestra$esp3, B = i))
## [1] 7.873408 3.515072 4.734132 3.484335 3.973367 3.782302 4.019471 4.132856
## [9] 4.211418

Ejemplos

Componentes principales: calificaciones en exámenes

Los datos marks (Mardia, Kent y Bibby, 1979) contienen los puntajes de 88 estudiantes en 5 pruebas: mecánica, vectores, álgebra, análisis y estadística. Cada renglón corresponde a la calificación de un estudiante en cada prueba.

marks <- read.csv("datos/marks.csv")
head(marks)
##   MECH VECT ALG ANL STAT
## 1   77   82  67  67   81
## 2   63   78  80  70   81
## 3   75   73  71  66   81
## 4   55   72  63  70   68
## 5   63   63  65  70   63
## 6   53   61  72  64   73

Entonces un análisis de componentes principales proseguiría como sigue:

pc_marks <- princomp(marks)
summary(pc_marks)
## Importance of components:
##                           Comp.1     Comp.2      Comp.3     Comp.4
## Standard deviation     26.061142 14.1355705 10.12760414 9.14706148
## Proportion of Variance  0.619115  0.1821424  0.09349705 0.07626893
## Cumulative Proportion   0.619115  0.8012575  0.89475453 0.97102347
##                            Comp.5
## Standard deviation     5.63807655
## Proportion of Variance 0.02897653
## Cumulative Proportion  1.00000000
loadings(pc_marks)
## 
## Loadings:
##      Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## MECH -0.505  0.749  0.300  0.296       
## VECT -0.368  0.207 -0.416 -0.783 -0.189
## ALG  -0.346        -0.145         0.924
## ANL  -0.451 -0.301 -0.597  0.518 -0.286
## STAT -0.535 -0.548  0.600 -0.176 -0.151
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings       1.0    1.0    1.0    1.0    1.0
## Proportion Var    0.2    0.2    0.2    0.2    0.2
## Cumulative Var    0.2    0.4    0.6    0.8    1.0
plot(pc_marks, type = "lines")

biplot(pc_marks)

Los cálculos de un análisis de componentes principales involucran la matriz de covarianzas empírica \(G\) (estimaciones plug-in)

\[G_{jk} = \frac{1}{88}\sum_{i=1}^88(x_{ij}-\bar{x_j})(x_{ik}-\bar{x_k})\]

para \(j,k=1,2,3,4,5\), y donde \(\bar{x_j} = \sum_{i=1}^88 x_{ij} / 88\) (la media de la i-ésima columna).

G <- cov(marks) * 87 / 88
G
##          MECH      VECT       ALG       ANL      STAT
## MECH 302.2934 125.77686 100.42510 105.06508 116.07076
## VECT 125.7769 170.87810  84.18957  93.59711  97.88688
## ALG  100.4251  84.18957 111.60318 110.83936 120.48567
## ANL  105.0651  93.59711 110.83936 217.87603 153.76808
## STAT 116.0708  97.88688 120.48567 153.76808 294.37177

Los pesos y las componentes principales no son mas que los eigenvalores y eigenvectores de la matriz de covarianzas \(G\), estos se calculan a través de una serie de de manipulaciones algebraicas que requieren cálculos del orden de p^3 (cuando G es una matriz de tamaño p\(\times\)p).

eigen_G <- eigen(G)
lambda <- eigen_G$values
v <- eigen_G$vectors
lambda
## [1] 679.18311 199.81435 102.56837  83.66873  31.78791
v
##            [,1]        [,2]       [,3]         [,4]        [,5]
## [1,] -0.5054457  0.74874751  0.2997888  0.296184264 -0.07939388
## [2,] -0.3683486  0.20740314 -0.4155900 -0.782888173 -0.18887639
## [3,] -0.3456612 -0.07590813 -0.1453182 -0.003236339  0.92392015
## [4,] -0.4511226 -0.30088849 -0.5966265  0.518139724 -0.28552169
## [5,] -0.5346501 -0.54778205  0.6002758 -0.175732020 -0.15123239
  1. Proponemos el siguiente modelo simple para puntajes correlacionados:

\[\textbf{x}_i = Q_i \textbf{v}\]

donde \(\textbf{x}_i\) es la tupla de calificaciones del i-ésimo estudiante, \(Q_i\) es un número que representa la habilidad del estudiante y \(\textbf{v}\) es un vector fijo con 5 números que aplica a todos los estudiantes. Si este modelo simple fuera cierto, entonces únicamente el \(\hat{\lambda}_1\) sería positivo y \(\textbf{v} = \hat{v}_1\). Sea \[\hat{\theta}=\sum_{i=1}^5\hat{\lambda}_i\] el modelo propuesto es equivalente a \(\hat{\theta}=1\), inculso si el modelo es correcto, no esperamos que \(\hat{\theta}\) sea exactamente uno pues hay ruido en los datos.

theta_hat <- lambda[1]/sum(lambda)
theta_hat
## [1] 0.619115

El valor de \(\hat{\theta}\) mide el porcentaje de la varianza explicada en la primer componente principal, ¿qué tan preciso es \(\hat{\theta}\)? La complejidad matemática en el cálculo de \(\hat{\theta}\) es irrelevante siempre y cuando podamos calcular \(\hat{\theta}^*\) para una muestra bootstrap, en esta caso una muestra bootsrtap es una base de datos de 88$$5 \(\textbf{X}^*\), donde las filas \(\textbf{x_i}^*\) de \(\textbf{X}^*\) son una muestra aleatoria de tamaño 88 de la verdadera matriz de datos.

pcBoot <- function(){
  muestra_boot <- sample_n(marks, size = 88, replace = TRUE)
  G <- cov(muestra_boot) * 87 / 88 
  eigen_G <- eigen(G)
  theta_hat <- eigen_G$values[1] / sum(eigen_G$values)
}
B <- 500
thetas_boot <- rdply(B, pcBoot)

Veamos un histograma de las replicaciones de \(\hat{\theta}\)

ggplot(thetas_boot, aes(x = V1)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.015, fill = "darkgray") + 
  geom_vline(aes(xintercept = mean(V1), color = "red"))

Estas tienen un error estándar

theta_se <- sd(thetas_boot$V1)
theta_se
## [1] 0.04639172

y media

mean(thetas_boot$V1)
## [1] 0.6156372

la media de las replicaciones es muy similar a la estimación \(\hat{\theta}\), esto indica que \(\hat{\theta}\) es cercano a insesgado. El intervalo de confianza estándar para el verdadero valor \({\theta}\) es:

\[\theta \in (\hat{\theta}- z^{(1-\alpha)}\cdot \hat{se}, \hat{\theta}-z^{(\alpha)}\cdot \hat{se})\] con un nivel de confianza de \(1-2\alpha\). Entonces:

# 0.68
theta_hat - theta_se
## [1] 0.5727233
theta_hat + theta_se
## [1] 0.6655068
# 0.9
theta_hat - qnorm(0.90) * theta_se
## [1] 0.5596617
theta_hat + qnorm(0.90) * theta_se
## [1] 0.6785684
  1. El eigenvetor \(\hat{v}_1\) correspondiente al mayor eigenvalor se conoce como primera componente de \(G\), supongamos que deseamos resumir la calificación de los estudiantes mediante un único número, entonces la mejor combinación lineal de los puntajes es

\[y_i = \sum_{k = 1}^5 \hat{v}_{1k}x_{ik}\]

esto es, la combinación lineal que utiliza las componentes de \(\hat{v}_1\) como ponderadores. Si queremos un resumen compuesto por dos números \((y_i,z_i)\), la segunda combinación lineal debería ser:

\[z_i = \sum_{k = 1}^5 \hat{v}_{2k}x_{ik}\]

Las componentes principales \(\hat{v}_1\) y \(\hat{v}_2\) son estadísticos, usa bootstrap para dar una medición de su variabilidad.