Ejercicio Bootstrap

library(tidyverse)
library(rsample)
library(GGally)
library(robust)
library(ggridges)
library(ggthemes)
library(magrittr)

El objetivo del ejercicio es comparar el ajuste lineal clásico (por mínimos cuadrados) con el ajuste robusto, en presencia de contaminación de la muestra. Evaluaremos la performance de los estimadores a través del bootstrap.

Modelo simulado.

Primero fijemos la semilla

set.seed(141414)

Para \(n=100\), simule el siguiente modelo.

Modelo que sigue el 90% de los datos

\[ Y=4+1.5X_{1}+8X_{2}-2X_{3}+\varepsilon \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1) \]

Modelo que sigue el 10% de los datos (los contaminados)

\[ Y=4+1.5X_{1}+8X_{2}-60X_{3}+\varepsilon \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2) \]

con

\(X_{1}\sim U[0,10],\)
\(X_{2}\sim U[0,10],\)
\(X_{3}\sim U[0,10],\)
\(\varepsilon\sim N (0,1),\)

independientes entre sí

Observemos que la única diferencia está en el coeficiente de la variable \(X_3\). El modelo (1) es el modelo verdadero.

n <- 100
x1 <- runif(n,0,10)
x2 <- runif(n,0,10)
x3 <- runif(n,0,10)
err <- rnorm(n,mean = 0, sd = 1)
# Con los valores de las X y el error calculados, construimos los dos modelos para y
y1 <- 4+1.5*x1[1:90]+8*x2[1:90]-2*x3[1:90]+err[1:90]
y2 <- 4+1.5*x1[91:100]+8*x2[91:100]-60*x3[91:100]+err[91:100]
# Jutnamos ambas y-es
y <- c(y1,y2)
# Va a resultar conveniente tener todo junto en un DF (el error no lo vamos a precisar)
df <- data.frame(x1,x2,x3,y)
#Podemos randomizar el orden de las filas, para que no queden todos los casos patológicos al final
df <- df[sample(nrow(df)),]
df

una buena forma de visualizar los datos sería con un diagrama de dispersión de Y respecto a las X

df %>% 
  gather(var,val,1:3) %>% 
  ggplot(.,aes(val,y,color=var))+
  geom_point()+
  facet_wrap(~var, scales = "free")+
  theme(legend.position = "none")

df %>% 
  ggplot(.,aes(y=y))+
  geom_point(aes(x3))+
  geom_abline(slope = -2,intercept = 4, color="green")+
  geom_abline(slope = -60,intercept = 4, color = "firebrick")

claramente se ven 10 casos con un comportamiento patológico.

En el segundo gráfico, podemos ver el detalle del comportamiento de las Y respecto a \(x_3\):

Notemos que las observaciones quedan por encima de las rectas trazadas por el efecto aditivo de \(x_1\) y \(x_2\)
Si no estuviera dicho efecto, las observaciones quedarían centradas sobre la recta correspondiente, con la dispersión agregada por el término del error.

Ajuste modelo lineal

Modelo clásico

modelo <- lm(y~x1+x2+x3,data = df)
summary(modelo)


Call:
lm(formula = y ~ x1 + x2 + x3, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-494.34    6.79   28.22   48.86   74.62 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   13.138     35.012   0.375   0.7083  
x1             2.746      3.931   0.699   0.4865  
x2             5.338      3.923   1.361   0.1768  
x3            -8.783      3.549  -2.475   0.0151 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 109.8 on 96 degrees of freedom
Multiple R-squared:  0.07943,   Adjusted R-squared:  0.05066 
F-statistic: 2.761 on 3 and 96 DF,  p-value: 0.0463

tidy(modelo,conf.int = TRUE,conf.level = 0.95)

glance(modelo)

au <- augment(modelo,df)

¿Son significativas las variables?

Modelo robusto

modelo_robusto <- lmRob(y~x1+x2+x3,data = df)

Max iteration for refinement reached.

summary(modelo_robusto)


Call:
lmRob(formula = y ~ x1 + x2 + x3, data = df)

Residuals:
      Min        1Q    Median        3Q       Max 
-563.6278   -0.7888   -0.1990    0.4478    3.5953 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.91814    0.35882   10.92   <2e-16 ***
x1           1.49676    0.03931   38.08   <2e-16 ***
x2           8.00335    0.03882  206.16   <2e-16 ***
x3          -1.98830    0.03478  -57.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9881 on 96 degrees of freedom
Multiple R-Squared: 0.7365 

Test for Bias:
            statistic p-value
M-estimate     4.7306  0.3161
LS-estimate    0.4702  0.9763

tidy(modelo_robusto)

glance(modelo_robusto)

au_rob <- augment(modelo_robusto,df)

En el modelo robusto ¿Son significativas las variables?

nota: Para el modelo lmRob no esta implementado confint(), por lo que no podemos calcular los intervalos de confianza en la función tidy()

Gráficos de los residuos

Notemos las diferencias en los gráficos de los residuos de los dos modelos:

ggplot(au, aes(.fitted, .resid)) +
  geom_point()+
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE)+
  labs(title= "OLS")

ggplot(au_rob, aes(.fitted, .resid)) +
  geom_point()+
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE)+
  labs(title= "Modelo Robusto")

para los valores estimados vs los residuos:

En ambos casos las observaciones patológicas toman valores extremos
lm(): Hay una estructura en los residuos de las observaciones no patológicas, hay una pendiente negativa entre los residuos y los valores estimados.
lmRob(): No hay pendiente, están todos cercanos a 0 (ajustan mejor!)

ggplot(au, aes(sample= .resid))+
  stat_qq()+
  geom_abline()+
  labs(title="QQ-plot OLS",
       x = "Distribución teórica",
       y= "Residuos")

ggplot(au_rob, aes(sample= .resid))+
  stat_qq()+
  geom_abline()+
  labs(title="QQ-plot Modelo Robusto",
       x = "Distribución teórica",
       y= "Residuos")

en los QQ-plots:

Nuevamente hay 10 observaciones patológicas que no siguen la distribución teórica en ningún caso
lm(): El resto de las observaciones tampoco sigue la distribución teórica, sino que en los cuantiles más altos hay una sobre-estimación de los residuos.
lmRob(): Hace un buen trabajo para las observaciones no patológicas

Bootstrap

Bootstrap¹

Extraemos una muestra de tamaño 100 con reemplazo de la muestra \({(X_{1i} , X_{2i} , X_{3i} , Y_i )}_{1≤i≤n}\) .
Ajustamos ambos estimadores (clásico y robusto) a cada muestra bootstrapeada.

Bootstrap casero

Primero hagamos un sampleo manual, para comprender los pasos.

Queremos tomar nuestro dataframe y construir 100 dataframes de la misma forma, donde cada observación se haya seleccionado de forma aleatoria de nuestro dataset original

set.seed(1234)
muestras_bootstrapeadas <- list()
#quiero construir 100 muestras bootstrpeadas. Podrían ser más
for (i in 1:100) { 
  #quiero que los dataset que muestree tengan el mismo tamaño que el dataset original
  muestra_i <- data.frame()
  index_i <- c()
  for (j in 1:nrow(df)) { 
    #eligo el número indice de la muestra de forma aleatoria entre los indices de mi set original
    index <- sample(rownames(df),1)
    #selecciono de mi df el elemento con el indice correspondiente
    muestra_unitaria <- df[index,]
    #agrego la observación a la muestra i
    muestra_i <- bind_rows(muestra_i,muestra_unitaria)
    #guardo el indice
    index_i <- c(index_i,index)
  }
  bootstrap_i <- list("df"=muestra_i, "indices" = index_i)
  muestras_bootstrapeadas[[i]] <- bootstrap_i
}

El resultado es una lista de 100 elementos, cada uno de los cuales es una muestra reconstruida a partir de la muestra original, donde tomamos los elementos con reposición: esto significa que podemos tomar más de una vez una observación.

Podemos ver esto revisando los índices de la primera muestra:

data.frame(table(muestras_bootstrapeadas[[1]]$indices)) %>% 
  arrange(-Freq)

por ejemplo, la observación 85 se tomo seis veces, la 56 cinco veces, etc.

Tidyverse al rescate

Sin embargo, hay formas más tidy de hacer un bootstrap. Para esto utilizamos la librería library(rsample)

muestras_bootstrapeadas <- bootstraps(df,times = 100)
muestras_bootstrapeadas

# Bootstrap sampling

Noten que el resultado es un df con dos columnas: un id y los splits que contienen para cada bootstrapeo un df.

Tomemos un elemento para ver qué tiene dentro

muestras_bootstrapeadas %>% 
  filter(id=="Bootstrap001") %$%
  splits[[1]][[1]]

data.frame(table(muestras_bootstrapeadas %>% 
  filter(id=="Bootstrap001") %$%
  splits[[1]][[2]])) %>% 
  arrange(-Freq)

Cada split contiene un dataframe con la muestra bootstrapeada, y un vector con los números de fila originales.

Si realizamos un table() podemos ver que algunos elementos tienen frecuencia mayor a 1. Esto es por realizamos un muestreo con reposición y por lo tanto hay una probabilidad de elegir más de una vez la misma muestra.

Este proceso replica el muestreo poblacional pero ahora consideramos a la muestra original como la población, de la cual tomamos n muestras.

Nos va a resultar conveniente envolver las funciones de forma tal que sólo necesiten un parámetro de entrada, para mapear

ajuste_lineal_simple <- function(df){
  lm(y~x1+x2+x3,data = df)
}
ajuste_lineal_robusto <- function(df){
  lmRob(y~x1+x2+x3,data = df)
}

Con esto podemos calcular el ajuste lineal y el ajuste robusto para cada una de las muestras.

muestras_bootstrapeadas <- muestras_bootstrapeadas %>% 
mutate(lm_simple = map(splits, ajuste_lineal_simple),
       lm_rob = map(splits, ajuste_lineal_robusto))
muestras_bootstrapeadas

# Bootstrap sampling

Para recuperar los parámetros de cada una de las estimaciones, vamos a realizar un map() sobre la función tidy()

parametros <- muestras_bootstrapeadas %>%
  gather(3:4) %>% 
  mutate(tdy = map(statistic,tidy)) %>% 
  unnest(tdy, .drop=TRUE)
parametros

Con todas estas estimaciones podemos graficar los resultados

parametros_verdaderos <- data_frame(term=c("(Intercept)","x1","x2","x3"),valor=c(4,1.5,8,-2))
ggplot(parametros, aes(estimate, fill = model))+
  geom_histogram(position= "dodge")+
  geom_vline(data=parametros_verdaderos,aes(xintercept = valor,color = "Verdadero"), size=1)+
  scale_color_manual(values = "darkolivegreen","")+
  theme_minimal()+
  scale_fill_gdocs()+
  theme(legend.position = "bottom")+
  facet_wrap(~term,scales = "free")

ggplot(parametros, aes(estimate,y=model, fill = model))+
  geom_density_ridges(alpha=.6)+
  theme_minimal()+
  geom_vline(data=parametros_verdaderos,aes(xintercept = valor,color = "Verdadero"), size=1)+
  scale_color_manual(values = "darkolivegreen","")+
  scale_fill_gdocs()+
  theme(legend.position = "bottom")+
  facet_wrap(~term,scales = "free")

Algunas conclusiones:

La distribución de los estimadores en el modelo robusto esta centrada en los parámetros verdaderos
La distribución de los estimadores en el modelo robusto esta Tiene mucha menos variabilidad
Particularmente para la variable patológica (\(x_3\)), la distribución del modelo lineal simple es especialmente mala, con el valor verdadero muy alejado del centro de masa.

diagrama basado en http://www.texample.net/tikz/examples/bootstrap-resampling/↩

LS0tCnRpdGxlOiAiRWplcmNpY2lvIEJvb3RzdHJhcCIKb3V0cHV0OiAKICBodG1sX25vdGVib29rOiAKICAgIGZpZ19oZWlnaHQ6IDEyCiAgICBmaWdfd2lkdGg6IDEwCiAgICB0aGVtZTogc3BhY2VsYWIKICAgIHRvYzogeWVzCiAgICB0b2NfZmxvYXQ6IHllcwotLS0KCgpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShyc2FtcGxlKQpsaWJyYXJ5KEdHYWxseSkKbGlicmFyeShyb2J1c3QpCmxpYnJhcnkoZ2dyaWRnZXMpCmxpYnJhcnkoZ2d0aGVtZXMpCmxpYnJhcnkobWFncml0dHIpCmBgYAoKRWwgb2JqZXRpdm8gZGVsIGVqZXJjaWNpbyBlcyBjb21wYXJhciBlbCBhanVzdGUgbGluZWFsIGNsw6FzaWNvIChwb3IgbcOtbmltb3MgY3VhZHJhZG9zKSBjb24gZWwgYWp1c3RlIHJvYnVzdG8sIGVuIHByZXNlbmNpYSBkZSBjb250YW1pbmFjacOzbiBkZSBsYSBtdWVzdHJhLiBFdmFsdWFyZW1vcyBsYSBwZXJmb3JtYW5jZSBkZSBsb3MgZXN0aW1hZG9yZXMgYSB0cmF2w6lzIGRlbCBib290c3RyYXAuCgojIyBNb2RlbG8gc2ltdWxhZG8uIAoKUHJpbWVybyBmaWplbW9zIGxhIHNlbWlsbGEKCmBgYHtyfQpzZXQuc2VlZCgxNDE0MTQpCmBgYAoKClBhcmEgJG49MTAwJCwgc2ltdWxlIGVsIHNpZ3VpZW50ZSBtb2RlbG8uCgoKX19Nb2RlbG8gcXVlIHNpZ3VlIGVsIDkwXCVfXyBkZSBsb3MgZGF0b3MKCiQkClk9NCsxLjVYX3sxfSs4WF97Mn0tMlhfezN9K1x2YXJlcHNpbG9uIApcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCAoMSkKJCQKCl9fTW9kZWxvIHF1ZSBzaWd1ZSBlbCAxMFwlX18gZGUgbG9zIGRhdG9zIChsb3MgY29udGFtaW5hZG9zKQoKJCQKWT00KzEuNVhfezF9KzhYX3syfS02MFhfezN9K1x2YXJlcHNpbG9uClwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcIFwgXCBcICgyKQokJAoKCmNvbgoKLSAkWF97MX1cc2ltIFVbMCwxMF0sJAotICRYX3syfVxzaW0gVVswLDEwXSwkCi0gJFhfezN9XHNpbSBVWzAsMTBdLCQKLSAkXHZhcmVwc2lsb25cc2ltIE4gKDAsMSksJAoKCmluZGVwZW5kaWVudGVzIGVudHJlIHPDrQoKT2JzZXJ2ZW1vcyBxdWUgbGEgw7puaWNhIGRpZmVyZW5jaWEgZXN0w6EgZW4gZWwgY29lZmljaWVudGUgZGUgbGEgdmFyaWFibGUgJFhfMyQuIF9FbCBtb2RlbG8gKDEpIGVzIGVsIG1vZGVsbyB2ZXJkYWRlcm8uXwoKCmBgYHtyfQpuIDwtIDEwMAp4MSA8LSBydW5pZihuLDAsMTApCngyIDwtIHJ1bmlmKG4sMCwxMCkKeDMgPC0gcnVuaWYobiwwLDEwKQplcnIgPC0gcm5vcm0obixtZWFuID0gMCwgc2QgPSAxKQojIENvbiBsb3MgdmFsb3JlcyBkZSBsYXMgWCB5IGVsIGVycm9yIGNhbGN1bGFkb3MsIGNvbnN0cnVpbW9zIGxvcyBkb3MgbW9kZWxvcyBwYXJhIHkKeTEgPC0gNCsxLjUqeDFbMTo5MF0rOCp4MlsxOjkwXS0yKngzWzE6OTBdK2VyclsxOjkwXQp5MiA8LSA0KzEuNSp4MVs5MToxMDBdKzgqeDJbOTE6MTAwXS02MCp4M1s5MToxMDBdK2Vycls5MToxMDBdCiMgSnV0bmFtb3MgYW1iYXMgeS1lcwp5IDwtIGMoeTEseTIpCiMgVmEgYSByZXN1bHRhciBjb252ZW5pZW50ZSB0ZW5lciB0b2RvIGp1bnRvIGVuIHVuIERGIChlbCBlcnJvciBubyBsbyB2YW1vcyBhIHByZWNpc2FyKQpkZiA8LSBkYXRhLmZyYW1lKHgxLHgyLHgzLHkpCgojUG9kZW1vcyByYW5kb21pemFyIGVsIG9yZGVuIGRlIGxhcyBmaWxhcywgcGFyYSBxdWUgbm8gcXVlZGVuIHRvZG9zIGxvcyBjYXNvcyBwYXRvbMOzZ2ljb3MgYWwgZmluYWwKZGYgPC0gZGZbc2FtcGxlKG5yb3coZGYpKSxdCmRmCmBgYAoKCnVuYSBidWVuYSBmb3JtYSBkZSB2aXN1YWxpemFyIGxvcyBkYXRvcyBzZXLDrWEgY29uIHVuIGRpYWdyYW1hIGRlIGRpc3BlcnNpw7NuIGRlIFkgcmVzcGVjdG8gYSBsYXMgWAoKYGBge3J9CmRmICU+JSAKICBnYXRoZXIodmFyLHZhbCwxOjMpICU+JSAKICBnZ3Bsb3QoLixhZXModmFsLHksY29sb3I9dmFyKSkrCiAgZ2VvbV9wb2ludCgpKwogIGZhY2V0X3dyYXAofnZhciwgc2NhbGVzID0gImZyZWUiKSsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpCgpkZiAlPiUgCiAgZ2dwbG90KC4sYWVzKHk9eSkpKwogIGdlb21fcG9pbnQoYWVzKHgzKSkrCiAgZ2VvbV9hYmxpbmUoc2xvcGUgPSAtMixpbnRlcmNlcHQgPSA0LCBjb2xvcj0iZ3JlZW4iKSsKICBnZW9tX2FibGluZShzbG9wZSA9IC02MCxpbnRlcmNlcHQgPSA0LCBjb2xvciA9ICJmaXJlYnJpY2siKQoKYGBgCgpjbGFyYW1lbnRlIHNlIHZlbiAxMCBjYXNvcyBjb24gdW4gY29tcG9ydGFtaWVudG8gcGF0b2zDs2dpY28uIAoKRW4gZWwgc2VndW5kbyBncsOhZmljbywgcG9kZW1vcyB2ZXIgZWwgZGV0YWxsZSBkZWwgY29tcG9ydGFtaWVudG8gZGUgbGFzIFkgcmVzcGVjdG8gYSAkeF8zJDoKCi0gTm90ZW1vcyBxdWUgbGFzIG9ic2VydmFjaW9uZXMgcXVlZGFuIHBvciBlbmNpbWEgZGUgbGFzIHJlY3RhcyB0cmF6YWRhcyBwb3IgZWwgZWZlY3RvIGFkaXRpdm8gZGUgJHhfMSQgeSAkeF8yJAotIFNpIG5vIGVzdHV2aWVyYSBkaWNobyBlZmVjdG8sIGxhcyBvYnNlcnZhY2lvbmVzIHF1ZWRhcsOtYW4gY2VudHJhZGFzIHNvYnJlIGxhIHJlY3RhIGNvcnJlc3BvbmRpZW50ZSwgY29uIGxhIGRpc3BlcnNpw7NuIGFncmVnYWRhIHBvciBlbCB0w6lybWlubyBkZWwgZXJyb3IuCgoKIyMgQWp1c3RlIG1vZGVsbyBsaW5lYWwKCiMjIyBNb2RlbG8gY2zDoXNpY28KCmBgYHtyfQoKbW9kZWxvIDwtIGxtKHl+eDEreDIreDMsZGF0YSA9IGRmKQoKc3VtbWFyeShtb2RlbG8pCnRpZHkobW9kZWxvLGNvbmYuaW50ID0gVFJVRSxjb25mLmxldmVsID0gMC45NSkKZ2xhbmNlKG1vZGVsbykKYXUgPC0gYXVnbWVudChtb2RlbG8sZGYpCgpgYGAKCgrCv1NvbiBzaWduaWZpY2F0aXZhcyBsYXMgdmFyaWFibGVzPwoKIyMjIE1vZGVsbyByb2J1c3RvCgpgYGB7cn0KbW9kZWxvX3JvYnVzdG8gPC0gbG1Sb2IoeX54MSt4Mit4MyxkYXRhID0gZGYpCgpzdW1tYXJ5KG1vZGVsb19yb2J1c3RvKQp0aWR5KG1vZGVsb19yb2J1c3RvKQpnbGFuY2UobW9kZWxvX3JvYnVzdG8pCgphdV9yb2IgPC0gYXVnbWVudChtb2RlbG9fcm9idXN0byxkZikKYGBgCgoKRW4gZWwgbW9kZWxvIHJvYnVzdG8gwr9Tb24gc2lnbmlmaWNhdGl2YXMgbGFzIHZhcmlhYmxlcz8KCl9ub3RhOiBQYXJhIGVsIG1vZGVsbyBsbVJvYiBubyBlc3RhIGltcGxlbWVudGFkbyBgY29uZmludCgpYCwgcG9yIGxvIHF1ZSBubyBwb2RlbW9zIGNhbGN1bGFyIGxvcyBpbnRlcnZhbG9zIGRlIGNvbmZpYW56YSBlbiBsYSBmdW5jacOzbiBgdGlkeSgpYF8gCgoKCiMjIyMgR3LDoWZpY29zIGRlIGxvcyByZXNpZHVvcwoKTm90ZW1vcyBsYXMgZGlmZXJlbmNpYXMgZW4gbG9zIGdyw6FmaWNvcyBkZSBsb3MgcmVzaWR1b3MgZGUgbG9zIGRvcyBtb2RlbG9zOgoKYGBge3J9CmdncGxvdChhdSwgYWVzKC5maXR0ZWQsIC5yZXNpZCkpICsKICBnZW9tX3BvaW50KCkrCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMCkgKwogIGdlb21fc21vb3RoKHNlID0gRkFMU0UpKwogIGxhYnModGl0bGU9ICJPTFMiKQoKCmdncGxvdChhdV9yb2IsIGFlcyguZml0dGVkLCAucmVzaWQpKSArCiAgZ2VvbV9wb2ludCgpKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDApICsKICBnZW9tX3Ntb290aChzZSA9IEZBTFNFKSsKICBsYWJzKHRpdGxlPSAiTW9kZWxvIFJvYnVzdG8iKQoKYGBgCgoKcGFyYSBsb3MgdmFsb3JlcyBlc3RpbWFkb3MgdnMgbG9zIHJlc2lkdW9zOgoKLSBFbiBhbWJvcyBjYXNvcyBsYXMgb2JzZXJ2YWNpb25lcyBwYXRvbMOzZ2ljYXMgdG9tYW4gdmFsb3JlcyBleHRyZW1vcwotIGBsbSgpYDogSGF5IHVuYSBlc3RydWN0dXJhIGVuIGxvcyByZXNpZHVvcyBkZSBsYXMgb2JzZXJ2YWNpb25lcyBubyBwYXRvbMOzZ2ljYXMsIGhheSB1bmEgcGVuZGllbnRlIG5lZ2F0aXZhIGVudHJlIGxvcyByZXNpZHVvcyB5IGxvcyB2YWxvcmVzIGVzdGltYWRvcy4KLSBgbG1Sb2IoKWA6IE5vIGhheSBwZW5kaWVudGUsIGVzdMOhbiB0b2RvcyBjZXJjYW5vcyBhIDAgKGFqdXN0YW4gbWVqb3IhKQoKCmBgYHtyfQoKZ2dwbG90KGF1LCBhZXMoc2FtcGxlPSAucmVzaWQpKSsKICBzdGF0X3FxKCkrCiAgZ2VvbV9hYmxpbmUoKSsKICBsYWJzKHRpdGxlPSJRUS1wbG90IE9MUyIsCiAgICAgICB4ID0gIkRpc3RyaWJ1Y2nDs24gdGXDs3JpY2EiLAogICAgICAgeT0gIlJlc2lkdW9zIikKCgpnZ3Bsb3QoYXVfcm9iLCBhZXMoc2FtcGxlPSAucmVzaWQpKSsKICBzdGF0X3FxKCkrCiAgZ2VvbV9hYmxpbmUoKSsKICBsYWJzKHRpdGxlPSJRUS1wbG90IE1vZGVsbyBSb2J1c3RvIiwKICAgICAgIHggPSAiRGlzdHJpYnVjacOzbiB0ZcOzcmljYSIsCiAgICAgICB5PSAiUmVzaWR1b3MiKQpgYGAKCgplbiBsb3MgUVEtcGxvdHM6CgotIE51ZXZhbWVudGUgaGF5IDEwIG9ic2VydmFjaW9uZXMgcGF0b2zDs2dpY2FzIHF1ZSBubyBzaWd1ZW4gbGEgZGlzdHJpYnVjacOzbiB0ZcOzcmljYSBlbiBuaW5nw7puIGNhc28KLSBgbG0oKWA6IEVsIHJlc3RvIGRlIGxhcyBvYnNlcnZhY2lvbmVzIHRhbXBvY28gc2lndWUgbGEgZGlzdHJpYnVjacOzbiB0ZcOzcmljYSwgc2lubyBxdWUgZW4gbG9zIGN1YW50aWxlcyBtw6FzIGFsdG9zIGhheSB1bmEgc29icmUtZXN0aW1hY2nDs24gZGUgbG9zIHJlc2lkdW9zLgotIGBsbVJvYigpYDogSGFjZSB1biBidWVuIHRyYWJham8gcGFyYSBsYXMgb2JzZXJ2YWNpb25lcyBubyBwYXRvbMOzZ2ljYXMKCgoKIyMgQm9vdHN0cmFwCgoKIVtCb290c3RyYXBeW2RpYWdyYW1hIGJhc2FkbyBlbiBodHRwOi8vd3d3LnRleGFtcGxlLm5ldC90aWt6L2V4YW1wbGVzL2Jvb3RzdHJhcC1yZXNhbXBsaW5nLwpdXShib290c3RyYXAtcmVzYW1wbGluZy5wbmcpCgoKLSBFeHRyYWVtb3MgdW5hIG11ZXN0cmEgZGUgdGFtYcOxbyAxMDAgY29uIHJlZW1wbGF6byBkZSBsYSBtdWVzdHJhICR7KFhfezFpfSAsIFhfezJpfSAsIFhfezNpfSAsIFlfaSApfV97MeKJpGniiaRufSQgLgoKLSBBanVzdGFtb3MgYW1ib3MgZXN0aW1hZG9yZXMgKGNsw6FzaWNvIHkgcm9idXN0bykgYSBjYWRhIG11ZXN0cmEgYm9vdHN0cmFwZWFkYS4KCgojIyMgQm9vdHN0cmFwIGNhc2VybwoKUHJpbWVybyBoYWdhbW9zIHVuIHNhbXBsZW8gbWFudWFsLCBwYXJhIGNvbXByZW5kZXIgbG9zIHBhc29zLgoKUXVlcmVtb3MgdG9tYXIgbnVlc3RybyBkYXRhZnJhbWUgeSBjb25zdHJ1aXIgMTAwIGRhdGFmcmFtZXMgZGUgbGEgbWlzbWEgZm9ybWEsIGRvbmRlIGNhZGEgb2JzZXJ2YWNpw7NuIHNlIGhheWEgc2VsZWNjaW9uYWRvIGRlIGZvcm1hIGFsZWF0b3JpYSBkZSBudWVzdHJvIGRhdGFzZXQgb3JpZ2luYWwKCmBgYHtyfQpzZXQuc2VlZCgxMjM0KQptdWVzdHJhc19ib290c3RyYXBlYWRhcyA8LSBsaXN0KCkKI3F1aWVybyBjb25zdHJ1aXIgMTAwIG11ZXN0cmFzIGJvb3RzdHJwZWFkYXMuIFBvZHLDrWFuIHNlciBtw6FzCmZvciAoaSBpbiAxOjEwMCkgeyAKICAjcXVpZXJvIHF1ZSBsb3MgZGF0YXNldCBxdWUgbXVlc3RyZWUgdGVuZ2FuIGVsIG1pc21vIHRhbWHDsW8gcXVlIGVsIGRhdGFzZXQgb3JpZ2luYWwKICBtdWVzdHJhX2kgPC0gZGF0YS5mcmFtZSgpCiAgaW5kZXhfaSA8LSBjKCkKICBmb3IgKGogaW4gMTpucm93KGRmKSkgeyAKICAgICNlbGlnbyBlbCBuw7ptZXJvIGluZGljZSBkZSBsYSBtdWVzdHJhIGRlIGZvcm1hIGFsZWF0b3JpYSBlbnRyZSBsb3MgaW5kaWNlcyBkZSBtaSBzZXQgb3JpZ2luYWwKICAgIGluZGV4IDwtIHNhbXBsZShyb3duYW1lcyhkZiksMSkKICAgICNzZWxlY2Npb25vIGRlIG1pIGRmIGVsIGVsZW1lbnRvIGNvbiBlbCBpbmRpY2UgY29ycmVzcG9uZGllbnRlCiAgICBtdWVzdHJhX3VuaXRhcmlhIDwtIGRmW2luZGV4LF0KICAgICNhZ3JlZ28gbGEgb2JzZXJ2YWNpw7NuIGEgbGEgbXVlc3RyYSBpCiAgICBtdWVzdHJhX2kgPC0gYmluZF9yb3dzKG11ZXN0cmFfaSxtdWVzdHJhX3VuaXRhcmlhKQogICAgI2d1YXJkbyBlbCBpbmRpY2UKICAgIGluZGV4X2kgPC0gYyhpbmRleF9pLGluZGV4KQogIH0KICBib290c3RyYXBfaSA8LSBsaXN0KCJkZiI9bXVlc3RyYV9pLCAiaW5kaWNlcyIgPSBpbmRleF9pKQogIG11ZXN0cmFzX2Jvb3RzdHJhcGVhZGFzW1tpXV0gPC0gYm9vdHN0cmFwX2kKfQoKYGBgCgoKRWwgcmVzdWx0YWRvIGVzIHVuYSBsaXN0YSBkZSAxMDAgZWxlbWVudG9zLCBjYWRhIHVubyBkZSBsb3MgY3VhbGVzIGVzIHVuYSBtdWVzdHJhIHJlY29uc3RydWlkYSBhIHBhcnRpciBkZSBsYSBtdWVzdHJhIG9yaWdpbmFsLCBkb25kZSB0b21hbW9zIGxvcyBlbGVtZW50b3MgX19jb24gcmVwb3NpY2nDs25fXzogZXN0byBzaWduaWZpY2EgcXVlIHBvZGVtb3MgdG9tYXIgbcOhcyBkZSB1bmEgdmV6IHVuYSBvYnNlcnZhY2nDs24uCgpQb2RlbW9zIHZlciBlc3RvIHJldmlzYW5kbyBsb3Mgw61uZGljZXMgZGUgbGEgcHJpbWVyYSBtdWVzdHJhOgoKCmBgYHtyIH0KZGF0YS5mcmFtZSh0YWJsZShtdWVzdHJhc19ib290c3RyYXBlYWRhc1tbMV1dJGluZGljZXMpKSAlPiUgCiAgYXJyYW5nZSgtRnJlcSkKYGBgCgpwb3IgZWplbXBsbywgbGEgb2JzZXJ2YWNpw7NuIDg1IHNlIHRvbW8gc2VpcyB2ZWNlcywgbGEgNTYgY2luY28gdmVjZXMsIGV0Yy4gCgojIyMgVGlkeXZlcnNlIGFsIHJlc2NhdGUKClNpbiBlbWJhcmdvLCBoYXkgZm9ybWFzIG3DoXMgX3RpZHlfIGRlIGhhY2VyIHVuIGJvb3RzdHJhcC4gUGFyYSBlc3RvIHV0aWxpemFtb3MgbGEgbGlicmVyw61hIGBsaWJyYXJ5KHJzYW1wbGUpYAoKYGBge3J9Cm11ZXN0cmFzX2Jvb3RzdHJhcGVhZGFzIDwtIGJvb3RzdHJhcHMoZGYsdGltZXMgPSAxMDApCm11ZXN0cmFzX2Jvb3RzdHJhcGVhZGFzCmBgYAoKTm90ZW4gcXVlIGVsIHJlc3VsdGFkbyBlcyB1biBkZiBjb24gZG9zIGNvbHVtbmFzOiB1biBpZCB5IGxvcyBfc3BsaXRzXyBxdWUgY29udGllbmVuIHBhcmEgY2FkYSBib290c3RyYXBlbyB1biBkZi4KClRvbWVtb3MgdW4gZWxlbWVudG8gcGFyYSB2ZXIgcXXDqSB0aWVuZSBkZW50cm8KCmBgYHtyfQoKbXVlc3RyYXNfYm9vdHN0cmFwZWFkYXMgJT4lIAogIGZpbHRlcihpZD09IkJvb3RzdHJhcDAwMSIpICUkJQogIHNwbGl0c1tbMV1dW1sxXV0KYGBgCgpgYGB7cn0KZGF0YS5mcmFtZSh0YWJsZShtdWVzdHJhc19ib290c3RyYXBlYWRhcyAlPiUgCiAgZmlsdGVyKGlkPT0iQm9vdHN0cmFwMDAxIikgJSQlCiAgc3BsaXRzW1sxXV1bWzJdXSkpICU+JSAKICBhcnJhbmdlKC1GcmVxKQoKYGBgCgoKQ2FkYSBzcGxpdCBjb250aWVuZSB1biBkYXRhZnJhbWUgY29uIGxhIG11ZXN0cmEgYm9vdHN0cmFwZWFkYSwgeSB1biB2ZWN0b3IgY29uIGxvcyBuw7ptZXJvcyBkZSBmaWxhIG9yaWdpbmFsZXMuIAoKU2kgcmVhbGl6YW1vcyB1biB0YWJsZSgpIHBvZGVtb3MgdmVyIHF1ZSBhbGd1bm9zIGVsZW1lbnRvcyB0aWVuZW4gZnJlY3VlbmNpYSBtYXlvciBhIDEuIEVzdG8gZXMgcG9yIHJlYWxpemFtb3MgdW4gbXVlc3RyZW8gX19jb24gcmVwb3NpY2nDs25fXyB5IHBvciBsbyB0YW50byBoYXkgdW5hIHByb2JhYmlsaWRhZCBkZSBlbGVnaXIgbcOhcyBkZSB1bmEgdmV6IGxhIG1pc21hIG11ZXN0cmEuIAoKRXN0ZSBwcm9jZXNvIF9fcmVwbGljYSBlbCBtdWVzdHJlbyBwb2JsYWNpb25hbF9fIHBlcm8gYWhvcmEgY29uc2lkZXJhbW9zIGEgX19sYSBtdWVzdHJhIG9yaWdpbmFsIGNvbW8gbGEgcG9ibGFjacOzbl9fLCBkZSBsYSBjdWFsIHRvbWFtb3MgbiBtdWVzdHJhcy4gCgoKTm9zIHZhIGEgcmVzdWx0YXIgY29udmVuaWVudGUgZW52b2x2ZXIgbGFzIGZ1bmNpb25lcyBkZSBmb3JtYSB0YWwgcXVlIHPDs2xvIG5lY2VzaXRlbiB1biBwYXLDoW1ldHJvIGRlIGVudHJhZGEsIHBhcmEgbWFwZWFyCmBgYHtyfQphanVzdGVfbGluZWFsX3NpbXBsZSA8LSBmdW5jdGlvbihkZil7CiAgbG0oeX54MSt4Mit4MyxkYXRhID0gZGYpCn0KCmFqdXN0ZV9saW5lYWxfcm9idXN0byA8LSBmdW5jdGlvbihkZil7CiAgbG1Sb2IoeX54MSt4Mit4MyxkYXRhID0gZGYpCn0KCmBgYAoKCkNvbiBlc3RvIHBvZGVtb3MgY2FsY3VsYXIgZWwgYWp1c3RlIGxpbmVhbCB5IGVsIGFqdXN0ZSByb2J1c3RvIHBhcmEgY2FkYSB1bmEgZGUgbGFzIG11ZXN0cmFzLiAKCmBgYHtyfQptdWVzdHJhc19ib290c3RyYXBlYWRhcyA8LSBtdWVzdHJhc19ib290c3RyYXBlYWRhcyAlPiUgCm11dGF0ZShsbV9zaW1wbGUgPSBtYXAoc3BsaXRzLCBhanVzdGVfbGluZWFsX3NpbXBsZSksCiAgICAgICBsbV9yb2IgPSBtYXAoc3BsaXRzLCBhanVzdGVfbGluZWFsX3JvYnVzdG8pKQoKCm11ZXN0cmFzX2Jvb3RzdHJhcGVhZGFzCmBgYAoKCgpQYXJhIHJlY3VwZXJhciBsb3MgcGFyw6FtZXRyb3MgZGUgY2FkYSB1bmEgZGUgbGFzIGVzdGltYWNpb25lcywgdmFtb3MgYSByZWFsaXphciB1biBgbWFwKClgIHNvYnJlIGxhIGZ1bmNpw7NuIGAgdGlkeSgpYAoKCmBgYHtyfQoKcGFyYW1ldHJvcyA8LSBtdWVzdHJhc19ib290c3RyYXBlYWRhcyAlPiUKICBnYXRoZXIoMzo0KSAlPiUgCiAgbXV0YXRlKHRkeSA9IG1hcChzdGF0aXN0aWMsdGlkeSkpICU+JSAKICB1bm5lc3QodGR5LCAuZHJvcD1UUlVFKQoKcGFyYW1ldHJvcwoKYGBgCgoKQ29uIHRvZGFzIGVzdGFzIGVzdGltYWNpb25lcyBwb2RlbW9zIGdyYWZpY2FyIGxvcyByZXN1bHRhZG9zCgpgYGB7cn0KcGFyYW1ldHJvc192ZXJkYWRlcm9zIDwtIGRhdGFfZnJhbWUodGVybT1jKCIoSW50ZXJjZXB0KSIsIngxIiwieDIiLCJ4MyIpLHZhbG9yPWMoNCwxLjUsOCwtMikpCgpnZ3Bsb3QocGFyYW1ldHJvcywgYWVzKGVzdGltYXRlLCBmaWxsID0gbW9kZWwpKSsKICBnZW9tX2hpc3RvZ3JhbShwb3NpdGlvbj0gImRvZGdlIikrCiAgZ2VvbV92bGluZShkYXRhPXBhcmFtZXRyb3NfdmVyZGFkZXJvcyxhZXMoeGludGVyY2VwdCA9IHZhbG9yLGNvbG9yID0gIlZlcmRhZGVybyIpLCBzaXplPTEpKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSAiZGFya29saXZlZ3JlZW4iLCIiKSsKICB0aGVtZV9taW5pbWFsKCkrCiAgc2NhbGVfZmlsbF9nZG9jcygpKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJib3R0b20iKSsKICBmYWNldF93cmFwKH50ZXJtLHNjYWxlcyA9ICJmcmVlIikKYGBgCgpgYGB7cn0KZ2dwbG90KHBhcmFtZXRyb3MsIGFlcyhlc3RpbWF0ZSx5PW1vZGVsLCBmaWxsID0gbW9kZWwpKSsKICBnZW9tX2RlbnNpdHlfcmlkZ2VzKGFscGhhPS42KSsKICB0aGVtZV9taW5pbWFsKCkrCiAgZ2VvbV92bGluZShkYXRhPXBhcmFtZXRyb3NfdmVyZGFkZXJvcyxhZXMoeGludGVyY2VwdCA9IHZhbG9yLGNvbG9yID0gIlZlcmRhZGVybyIpLCBzaXplPTEpKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSAiZGFya29saXZlZ3JlZW4iLCIiKSsKICBzY2FsZV9maWxsX2dkb2NzKCkrCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gImJvdHRvbSIpKwogIGZhY2V0X3dyYXAofnRlcm0sc2NhbGVzID0gImZyZWUiKQpgYGAKCkFsZ3VuYXMgY29uY2x1c2lvbmVzOgoKLSBMYSBkaXN0cmlidWNpw7NuIGRlIGxvcyBlc3RpbWFkb3JlcyBlbiBlbCBtb2RlbG8gcm9idXN0byBlc3RhIF9fY2VudHJhZGEgZW4gbG9zIHBhcsOhbWV0cm9zIHZlcmRhZGVyb3NfXwotIExhIGRpc3RyaWJ1Y2nDs24gZGUgbG9zIGVzdGltYWRvcmVzIGVuIGVsIG1vZGVsbyByb2J1c3RvIGVzdGEgX19UaWVuZSBtdWNoYSBtZW5vcyB2YXJpYWJpbGlkYWRfXwotIFBhcnRpY3VsYXJtZW50ZSBwYXJhIGxhIHZhcmlhYmxlIHBhdG9sw7NnaWNhICgkeF8zJCksIGxhIGRpc3RyaWJ1Y2nDs24gZGVsIG1vZGVsbyBsaW5lYWwgc2ltcGxlIGVzIGVzcGVjaWFsbWVudGUgbWFsYSwgY29uIGVsIHZhbG9yIHZlcmRhZGVybyBtdXkgYWxlamFkbyBkZWwgY2VudHJvIGRlIG1hc2EuCgoKCgo=