Objetivos

  1. Regresion Lineal
  2. Regresion Lineal multiple
  3. Otras consideraciones de modelos de regresion lineal

Introduccion

Se utilizara el dataset Advertising

data <- read.csv('Advertising.csv')
data <- data[,-1]
attach(data)
names(data)
[1] "TV"        "radio"     "newspaper" "sales"    
summary(data)
       TV             radio          newspaper          sales      
 Min.   :  0.70   Min.   : 0.000   Min.   :  0.30   Min.   : 1.60  
 1st Qu.: 74.38   1st Qu.: 9.975   1st Qu.: 12.75   1st Qu.:10.38  
 Median :149.75   Median :22.900   Median : 25.75   Median :12.90  
 Mean   :147.04   Mean   :23.264   Mean   : 30.55   Mean   :14.02  
 3rd Qu.:218.82   3rd Qu.:36.525   3rd Qu.: 45.10   3rd Qu.:17.40  
 Max.   :296.40   Max.   :49.600   Max.   :114.00   Max.   :27.00

Regresion lineal simple

Es de la forma \[Y \approx \beta_0 + \beta_1X + \epsilon_i\]

\(Y\) es la respuesta, o variable dependiente, y \(X\) es el predictor. En el caso puntual de la base de publicidad sobre la que estamos trabajando el modelo de regresion sería de la forma \[ sales \approx \beta_0 + \beta_1TV \] Donde \(\beta_0\) (ordenada al origen) y \(\beta_1\) (pendiente) son los coeficientes de regresión.

A continuación, generamos un diagrama de dispersión con los datos de Sales e inversión en publicidad por TV.

data %>% select(sales, TV) %>%
  ggplot(aes(x=TV, y=sales)) +
  geom_point()

data %>% select(sales, TV) %>%
  ggplot(aes(x=TV, y=sales)) +
  geom_point() + 
  stat_quantile(quantiles = seq(0.05, 0.95, by=0.05))

A priori, podemos obervar que existen infinitas rectas, o modelos, que pueden representar la información de nuestro dataset. Con MCO (minimos cuadrados ordinarios) minimizamos el error cuadrático entre nuestra aproximación \(\hat{Y}\) y los datos de la muestra \(\hat{Y}\). De esta manera, se encuentra la recta con coeficientes \(\hat{\beta_0}\) y \(\hat{\beta_1}\) que minimizan el error.

\[ min{RSS} = \sum{e_i^2} = \sum{(Y_i - \hat{Y_i})^2} = \sum{(Y_i - \hat{\beta_0} - \hat{\beta_1}X_i)^2}\]

rss <- lm(sales~TV, data=data)
summary(rss)

Call:
lm(formula = sales ~ TV, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.3860 -1.9545 -0.1913  2.0671  7.2124 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept) 7.032594   0.457843   15.36 <0.0000000000000002 ***
TV          0.047537   0.002691   17.67 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.259 on 198 degrees of freedom
Multiple R-squared:  0.6119,    Adjusted R-squared:  0.6099 
F-statistic: 312.1 on 1 and 198 DF,  p-value: < 0.00000000000000022

En \(R\), utilizamos la funcion \(lm()\) para hallar los coeficientes de la recta de regresión.

fit <- lm(sales~TV, data=data)
summary(fit)

Call:
lm(formula = sales ~ TV, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.3860 -1.9545 -0.1913  2.0671  7.2124 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept) 7.032594   0.457843   15.36 <0.0000000000000002 ***
TV          0.047537   0.002691   17.67 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.259 on 198 degrees of freedom
Multiple R-squared:  0.6119,    Adjusted R-squared:  0.6099 
F-statistic: 312.1 on 1 and 198 DF,  p-value: < 0.00000000000000022

En \(fit\$residuals\) encontramos los residuos de cada observación \(e_i = Y_i-\hat{Y_i}\)

fit$residuals[1:10]
         1          2          3          4          5          6          7 
 4.1292255  1.2520260  1.4497762  4.2656054 -2.7272181 -0.2461623  2.0340496 
         8          9         10 
 0.4535023 -2.6414087 -5.9304143

Diagrama de disperción partiendo de modelo

En este caso, podemos partir de un modelo establecido y, generando valores aleatorios con una determinada distribución de probabilida, sumar las variaciones al modelo para tener una diagrama de disperción y poder hacer el ejercicio inverso.

sim_data <- function(){
  x <- runif(30, min=-3, max=3)
  epsilon <- rnorm(n=length(x),m=0,sd=1)
  y <- -2+3*x+epsilon
  return(data.frame(x,y))
}

dataset1 <- sim_data()
dataset2 <- sim_data()
dataset3 <- sim_data()
dataset4 <- sim_data()

par(mfrow=c(2,2))
plot(dataset1$x, dataset1$y, main="dataset 1")
abline(lm(dataset1$y~dataset1$x, data=data), col="red")
plot(dataset1$x, dataset1$y, main="dataset 2")
abline(lm(dataset2$y~dataset2$x, data=data), col="red")
plot(dataset1$x, dataset1$y, main="dataset 3")
abline(lm(dataset3$y~dataset3$x, data=data), col="red")
plot(dataset1$x, dataset1$y, main="dataset 4")
abline(lm(dataset4$y~dataset4$x, data=data), col="red")

Calculamos coeficientes de regresión. Del modelo planteado, sabemos que deberian ser \(\beta_0 = -2\) y \(\beta_1 = 3\)

sim_regr_1 <- lm(dataset1$y~dataset1$x, data=data)
sim_regr_2 <- lm(dataset2$y~dataset2$x, data=data)
sim_regr_3 <- lm(dataset3$y~dataset3$x, data=data)
sim_regr_4 <- lm(dataset4$y~dataset4$x, data=data)
Dataset \(\beta_0\) \(\beta_1\)
Dataset1 -1.782959 2.951345
Dataset2 -1.9549081 2.7171619
Dataset3 -2.087374 2.9309085
Dataset4 -2.0500293 3.1399282

Vemos que, en todos los casos, tenemos un error asociado para los caluclos de los coeficientes de regresión.

Analizamos los errores estándar \(SE\) para el primer modelo.

  • Error estandar \(\beta_0\): 0.2337201
  • Error estandar \(\beta_1\): 0.1277985

Podemos establecer intervalos de congianza al \(95\%\) y \(99\%\) respectivamente para los coeficientes de regresión

\[ IC_{95\%} = [ \hat{\beta_i} - 2SE(\hat{\beta_i}), \hat{\beta_i} + 2SE(\hat{\beta_i}) ] \] \[ IC_{99\%} = [ \hat{\beta_i} - 3SE(\hat{\beta_i}), \hat{\beta_i} + 3SE(\hat{\beta_i}) ] \]

Parametro Intervalo
\(\beta_0\) [-2.2503992,-1.3155189]
\(\beta_1\) [2.6957479,3.206942]

Relacion entre X e Y en MCO

Se realiza un test de hipótesis nula para el coeficiente de regresión \(\beta_1\) que relaciona los dos parámetros del modelo. En este caso puntual serían \(sales\) y \(TV\). \(H_0\) es la hipótesis nula y \(H_a\) es la hipótesis alternativa.

  • \(H_0:\) No hay relacion entre \(sales\) y \(TV\)
  • \(H_a:\) Hay relacion entre \(sales\) y \(TV\)

Si \(\beta_1 = 0\) entonces no existe relacion entre ambos parámetros.

  • \(H_0: \beta_1 = 0\)
  • \(H_a: \beta_1 \neq 0\)

Si \(\beta_1\) es estadísticamente significativo entonces decimos que explica a la varianza de la viarable dependiente \(sales\).

Parametro P-valor
\(\beta_0\) 0
\(\beta_1\) 0

Regresion Lineal Multiple

Podemos analizar gráficamente la relación entre las distintas variables de la base.

plot(data)

rlm <- lm(sales ~ radio + TV + newspaper, data=data)
summary(rlm)

Call:
lm(formula = sales ~ radio + TV + newspaper, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8277 -0.8908  0.2418  1.1893  2.8292 

Coefficients:
             Estimate Std. Error t value            Pr(>|t|)    
(Intercept)  2.938889   0.311908   9.422 <0.0000000000000002 ***
radio        0.188530   0.008611  21.893 <0.0000000000000002 ***
TV           0.045765   0.001395  32.809 <0.0000000000000002 ***
newspaper   -0.001037   0.005871  -0.177                0.86    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.686 on 196 degrees of freedom
Multiple R-squared:  0.8972,    Adjusted R-squared:  0.8956 
F-statistic: 570.3 on 3 and 196 DF,  p-value: < 0.00000000000000022

\(\textbf{Conclusiones}\): Los coeficientes tienen la misma unidad que la variable independiente, ya sea esta \(\%\), \(\$\), u otra. Segun la regresión realizada, una inversión de 1 unidad en radio, es lo que generará más aumento en \(sales\). Tambien podemos que \(newspaper\) no es significativo estadísticamente en términos de la variación de \(sales\).

Realizamos la regresion lineal quitando a la variable \(newspaper\) por ser estadísticamente poco significativa. Se observa que los coeficientes no sufren variaciones importantes.

rlm2 <- lm(sales ~ radio + TV , data=data)
summary(rlm2)

Call:
lm(formula = sales ~ radio + TV, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.7977 -0.8752  0.2422  1.1708  2.8328 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept)  2.92110    0.29449   9.919 <0.0000000000000002 ***
radio        0.18799    0.00804  23.382 <0.0000000000000002 ***
TV           0.04575    0.00139  32.909 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.681 on 197 degrees of freedom
Multiple R-squared:  0.8972,    Adjusted R-squared:  0.8962 
F-statistic: 859.6 on 2 and 197 DF,  p-value: < 0.00000000000000022

Regresion Lineal con variables dicotómicas

credit <- read.csv('Credit.csv')

names(credit) <- tolower(names(credit))
credit <- credit[,-1]
head(credit)

Se modela el balance respecto de la variable dicotomica \(\textbf{gender}\). La categoría base será \(\textbf{female}\).

fit4 <- lm(balance ~ gender, data=credit)
summary(fit4)$coefficients
             Estimate Std. Error    t value
(Intercept) 529.53623   31.98819 16.5541153
genderMale  -19.73312   46.05121 -0.4285039
                                                           Pr(>|t|)
(Intercept) 0.00000000000000000000000000000000000000000000003312981
genderMale  0.66851610550271733934835083346115425229072570800781250

Podemos interpretar a \(\beta_0\) como el balance del submodelo que representa a las mujeres, y \(\beta_1\) como la varianza que tiene el balance de una persona por ser hombre.

\[ gender_i = \left\{\begin{matrix} 0 & Female\\ 1 & Male \end{matrix}\right.\] \[ balance = \beta_0 + \beta_1 \cdot gender = \left\{\begin{matrix} \beta_0 + \beta_1 & male\\ \beta_1 & female \end{matrix}\right. \]

Y si queremos cambiar la categoría base o los valores de las variables dicotómicas?

\[ gender_i = \left\{\begin{matrix} -1 & female\\ 1 & male \end{matrix}\right. \]

Esto deja una nueva expresión del modelo: \[ balance = \beta_0 + \beta_1 \cdot gender = \left\{\begin{matrix} \beta_0 + \beta_1 & male\\ \beta_0 - \beta_1 & female \end{matrix}\right. \]

credit$coded_gender <- ifelse(credit$gender == "Female",-1,1)
fit5 <- lm(balance~coded_gender, data=credit)
summary(fit5)$coefficients
               Estimate Std. Error    t value
(Intercept)  519.669670   23.02561 22.5692080
coded_gender  -9.866562   23.02561 -0.4285039
                                                                                      Pr(>|t|)
(Intercept)  0.0000000000000000000000000000000000000000000000000000000000000000000000003191322
coded_gender 0.6685161055027173393483508334611542522907257080078125000000000000000000000000000

En este caso, los valores de los coeficientes no cambian, solo su interpretación. Ya que \(\beta_1\) en este caso será varianza positica para género mujer y varianza negativa para género hombre.

Variables Multinomiales

Este tipo de variables discretas, que pueden tomar finitos valores, se dicotomizan.

Ejemplo: Tenemos una variable llamada \(\textit{ethniticy}\) con 3 categorías.

#Obtenemos valores que puede tomar la variable
unique(credit$ethnicity)
[1] Caucasian        Asian            African American
Levels: African American Asian Caucasian

\[ \beta_1 = asian_i = \left\{\begin{matrix} 0 & no\\ 1 & yes \end{matrix}\right.\]

\[ \beta_2 = Caucasian_i = \left\{\begin{matrix} 0 & no\\ 1 & yes \end{matrix}\right. \]

credit$caucasian <- ifelse(credit$ethnicity == "Caucasian",1,0)
credit$asian <- ifelse(credit$ethnicity == "Asian",1,0)

\[ balance = \beta_0 + \beta_1 \cdot asian + \beta_2 \cdot caucasian = \left\{\begin{matrix} \beta_0 + \beta_1 & asian\\ \beta_0 + \beta_2 & caucasian\\ \beta_0 & african\_american \end{matrix}\right. \]

fit6 <- lm(balance ~ asian + caucasian, data=credit)
summary(fit6)

Call:
lm(formula = balance ~ asian + caucasian, data = credit)

Residuals:
    Min      1Q  Median      3Q     Max 
-531.00 -457.08  -63.25  339.25 1480.50 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept)   531.00      46.32  11.464 <0.0000000000000002 ***
asian         -18.69      65.02  -0.287               0.774    
caucasian     -12.50      56.68  -0.221               0.826    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 460.9 on 397 degrees of freedom
Multiple R-squared:  0.0002188, Adjusted R-squared:  -0.004818 
F-statistic: 0.04344 on 2 and 397 DF,  p-value: 0.9575
LS0tCnRpdGxlOiAiUmVncmVzaW9uIExpbmVhbCIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCgpgYGB7ciwgZWNobz1GQUxTRSxldmFsPVRSVUUsd2FybmluZz1GQUxTRX0KI2luc3RhbGwucGFja2FnZXMoJ0lTTFInKQojaW5zdGFsbC5wYWNrYWdlcygnTUFTUycpCiNpbnN0YWxsLnBhY2thZ2VzKCdjYXInKQojaW5zdGFsbC5wYWNrYWdlcygncXVhbnRyZWcnKQojaW5zdGFsbC5wYWNrYWdlcygnc3BhcnNlTScpCgpsaWJyYXJ5KFNwYXJzZU0pCmxpYnJhcnkoY2FyKQpsaWJyYXJ5KE1BU1MpCmxpYnJhcnkoSVNMUikKbGlicmFyeShsbXRlc3QpCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkoZ2dwbG90MikKCm9wdGlvbnMoc2NpcGVuID0gOTk5OTkpICMgY2FtYmlvIGZvcm1hdG8gZGUgc2hvdyBkZWNpbWFscwpzZXR3ZCgnL2hvbWUvamdvbnphbGV6L2Rldi9kYXRhU2NpZW5jZScpCmBgYAoKIyMgT2JqZXRpdm9zCgoxLiBSZWdyZXNpb24gTGluZWFsCjIuIFJlZ3Jlc2lvbiBMaW5lYWwgbXVsdGlwbGUKMy4gT3RyYXMgY29uc2lkZXJhY2lvbmVzIGRlIG1vZGVsb3MgZGUgcmVncmVzaW9uIGxpbmVhbAoKCiMjIEludHJvZHVjY2lvbgoKU2UgdXRpbGl6YXJhIGVsIGRhdGFzZXQgYEFkdmVydGlzaW5nYAoKYGBge3IsZXZhbD1GQUxTRX0KZGF0YSA8LSByZWFkLmNzdignQWR2ZXJ0aXNpbmcuY3N2JykKZGF0YSA8LSBkYXRhWywtMV0KYXR0YWNoKGRhdGEpCmBgYAoKCmBgYHtyfQpuYW1lcyhkYXRhKQpgYGAKCmBgYHtyfQpzdW1tYXJ5KGRhdGEpCmBgYAoKCiMjIFJlZ3Jlc2lvbiBsaW5lYWwgc2ltcGxlCgpFcyBkZSBsYSBmb3JtYSAkJFkgXGFwcHJveCBcYmV0YV8wICsgXGJldGFfMVggKyBcZXBzaWxvbl9pJCQKCiRZJCBlcyBsYSByZXNwdWVzdGEsIG8gdmFyaWFibGUgZGVwZW5kaWVudGUsIHkgJFgkIGVzIGVsIHByZWRpY3Rvci4gRW4gZWwgY2FzbyBwdW50dWFsIGRlIGxhIGJhc2UgZGUgcHVibGljaWRhZCBzb2JyZSBsYSBxdWUgZXN0YW1vcyB0cmFiYWphbmRvIGVsIG1vZGVsbyBkZSByZWdyZXNpb24gc2Vyw61hIGRlIGxhIGZvcm1hICQkIHNhbGVzIFxhcHByb3ggXGJldGFfMCArIFxiZXRhXzFUViAkJApEb25kZSAkXGJldGFfMCQgKG9yZGVuYWRhIGFsIG9yaWdlbikgeSAkXGJldGFfMSQgKHBlbmRpZW50ZSkgc29uIGxvcyBjb2VmaWNpZW50ZXMgZGUgcmVncmVzacOzbi4KCkEgY29udGludWFjacOzbiwgZ2VuZXJhbW9zIHVuIGRpYWdyYW1hIGRlIGRpc3BlcnNpw7NuIGNvbiBsb3MgZGF0b3MgZGUgU2FsZXMgZSBpbnZlcnNpw7NuIGVuIHB1YmxpY2lkYWQgcG9yIFRWLgoKYGBge3J9CmRhdGEgJT4lIHNlbGVjdChzYWxlcywgVFYpICU+JQogIGdncGxvdChhZXMoeD1UViwgeT1zYWxlcykpICsKICBnZW9tX3BvaW50KCkKYGBgCgoKCmBgYHtyfQpkYXRhICU+JSBzZWxlY3Qoc2FsZXMsIFRWKSAlPiUKICBnZ3Bsb3QoYWVzKHg9VFYsIHk9c2FsZXMpKSArCiAgZ2VvbV9wb2ludCgpICsgCiAgc3RhdF9xdWFudGlsZShxdWFudGlsZXMgPSBzZXEoMC4wNSwgMC45NSwgYnk9MC4wNSkpCmBgYAoKQSBwcmlvcmksIHBvZGVtb3Mgb2JlcnZhciBxdWUgZXhpc3RlbiBpbmZpbml0YXMgcmVjdGFzLCBvIG1vZGVsb3MsIHF1ZSBwdWVkZW4gcmVwcmVzZW50YXIgbGEgaW5mb3JtYWNpw7NuIGRlIG51ZXN0cm8gZGF0YXNldC4gQ29uIE1DTyAobWluaW1vcyBjdWFkcmFkb3Mgb3JkaW5hcmlvcykgbWluaW1pemFtb3MgZWwgZXJyb3IgY3VhZHLDoXRpY28gZW50cmUgbnVlc3RyYSBhcHJveGltYWNpw7NuICRcaGF0e1l9JCB5IGxvcyBkYXRvcyBkZSBsYSBtdWVzdHJhICRcaGF0e1l9JC4gRGUgZXN0YSBtYW5lcmEsIHNlIGVuY3VlbnRyYSBsYSByZWN0YSBjb24gY29lZmljaWVudGVzICRcaGF0e1xiZXRhXzB9JCB5ICRcaGF0e1xiZXRhXzF9JCBxdWUgbWluaW1pemFuIGVsIGVycm9yLgoKJCQgbWlue1JTU30gPSBcc3Vte2VfaV4yfSA9IFxzdW17KFlfaSAtIFxoYXR7WV9pfSleMn0gPSBcc3VteyhZX2kgLSBcaGF0e1xiZXRhXzB9IC0gXGhhdHtcYmV0YV8xfVhfaSleMn0kJAoKCmBgYHtyfQpyc3MgPC0gbG0oc2FsZXN+VFYsIGRhdGE9ZGF0YSkKc3VtbWFyeShyc3MpCmBgYAoKRW4gJFIkLCB1dGlsaXphbW9zIGxhIGZ1bmNpb24gJGxtKCkkIHBhcmEgaGFsbGFyIGxvcyBjb2VmaWNpZW50ZXMgZGUgbGEgcmVjdGEgZGUgcmVncmVzacOzbi4KCmBgYHtyfQpmaXQgPC0gbG0oc2FsZXN+VFYsIGRhdGE9ZGF0YSkKc3VtbWFyeShmaXQpCmBgYAoKRW4gJGZpdFwkcmVzaWR1YWxzJCBlbmNvbnRyYW1vcyBsb3MgcmVzaWR1b3MgZGUgY2FkYSBvYnNlcnZhY2nDs24gJGVfaSA9IFlfaS1caGF0e1lfaX0kCmBgYHtyfQpmaXQkcmVzaWR1YWxzWzE6MTBdCmBgYAoKKiBUb3RhbCByZXNpZHVvczogYHIgbGVuZ3RoKGZpdCRyZXNpZHVhbHMpYAoqIFRvdGFsIG9ic2VydmFjaW9uZXM6IGByIG5yb3coZGF0YSlgCgoKIyMjIERpYWdyYW1hIGRlIGRpc3BlcmNpw7NuIHBhcnRpZW5kbyBkZSBtb2RlbG8KCkVuIGVzdGUgY2FzbywgcG9kZW1vcyBwYXJ0aXIgZGUgdW4gbW9kZWxvIGVzdGFibGVjaWRvIHksIGdlbmVyYW5kbyB2YWxvcmVzIGFsZWF0b3Jpb3MgY29uIHVuYSBkZXRlcm1pbmFkYSBkaXN0cmlidWNpw7NuIGRlIHByb2JhYmlsaWRhLCBzdW1hciBsYXMgdmFyaWFjaW9uZXMgYWwgbW9kZWxvIHBhcmEgdGVuZXIgdW5hIGRpYWdyYW1hIGRlIGRpc3BlcmNpw7NuIHkgcG9kZXIgaGFjZXIgZWwgZWplcmNpY2lvIGludmVyc28uCgpgYGB7cn0Kc2ltX2RhdGEgPC0gZnVuY3Rpb24oKXsKICB4IDwtIHJ1bmlmKDMwLCBtaW49LTMsIG1heD0zKQogIGVwc2lsb24gPC0gcm5vcm0obj1sZW5ndGgoeCksbT0wLHNkPTEpCiAgeSA8LSAtMiszKngrZXBzaWxvbgogIHJldHVybihkYXRhLmZyYW1lKHgseSkpCn0KCmRhdGFzZXQxIDwtIHNpbV9kYXRhKCkKZGF0YXNldDIgPC0gc2ltX2RhdGEoKQpkYXRhc2V0MyA8LSBzaW1fZGF0YSgpCmRhdGFzZXQ0IDwtIHNpbV9kYXRhKCkKCnBhcihtZnJvdz1jKDIsMikpCnBsb3QoZGF0YXNldDEkeCwgZGF0YXNldDEkeSwgbWFpbj0iZGF0YXNldCAxIikKYWJsaW5lKGxtKGRhdGFzZXQxJHl+ZGF0YXNldDEkeCwgZGF0YT1kYXRhKSwgY29sPSJyZWQiKQpwbG90KGRhdGFzZXQxJHgsIGRhdGFzZXQxJHksIG1haW49ImRhdGFzZXQgMiIpCmFibGluZShsbShkYXRhc2V0MiR5fmRhdGFzZXQyJHgsIGRhdGE9ZGF0YSksIGNvbD0icmVkIikKcGxvdChkYXRhc2V0MSR4LCBkYXRhc2V0MSR5LCBtYWluPSJkYXRhc2V0IDMiKQphYmxpbmUobG0oZGF0YXNldDMkeX5kYXRhc2V0MyR4LCBkYXRhPWRhdGEpLCBjb2w9InJlZCIpCnBsb3QoZGF0YXNldDEkeCwgZGF0YXNldDEkeSwgbWFpbj0iZGF0YXNldCA0IikKYWJsaW5lKGxtKGRhdGFzZXQ0JHl+ZGF0YXNldDQkeCwgZGF0YT1kYXRhKSwgY29sPSJyZWQiKQpgYGAKCkNhbGN1bGFtb3MgY29lZmljaWVudGVzIGRlIHJlZ3Jlc2nDs24uIERlbCBtb2RlbG8gcGxhbnRlYWRvLCBzYWJlbW9zIHF1ZSBkZWJlcmlhbiBzZXIgJFxiZXRhXzAgPSAtMiQgeSAkXGJldGFfMSA9IDMkCgpgYGB7cn0Kc2ltX3JlZ3JfMSA8LSBsbShkYXRhc2V0MSR5fmRhdGFzZXQxJHgsIGRhdGE9ZGF0YSkKc2ltX3JlZ3JfMiA8LSBsbShkYXRhc2V0MiR5fmRhdGFzZXQyJHgsIGRhdGE9ZGF0YSkKc2ltX3JlZ3JfMyA8LSBsbShkYXRhc2V0MyR5fmRhdGFzZXQzJHgsIGRhdGE9ZGF0YSkKc2ltX3JlZ3JfNCA8LSBsbShkYXRhc2V0NCR5fmRhdGFzZXQ0JHgsIGRhdGE9ZGF0YSkKYGBgCgoKRGF0YXNldCB8ICRcYmV0YV8wJCB8ICRcYmV0YV8xJCB8IAotLS0tLS0tLXwtLS0tLS0tLS0tLXwtLS0tLS0tLS0tLXwKRGF0YXNldDF8YHIgc2ltX3JlZ3JfMSRjb2VmZmljaWVudHNbMV1gfGByIHNpbV9yZWdyXzEkY29lZmZpY2llbnRzWzJdYHwKRGF0YXNldDJ8YHIgc2ltX3JlZ3JfMiRjb2VmZmljaWVudHNbMV1gfGByIHNpbV9yZWdyXzIkY29lZmZpY2llbnRzWzJdYHwKRGF0YXNldDN8YHIgc2ltX3JlZ3JfMyRjb2VmZmljaWVudHNbMV1gfGByIHNpbV9yZWdyXzMkY29lZmZpY2llbnRzWzJdYHwKRGF0YXNldDR8YHIgc2ltX3JlZ3JfNCRjb2VmZmljaWVudHNbMV1gfGByIHNpbV9yZWdyXzQkY29lZmZpY2llbnRzWzJdYHwKClZlbW9zIHF1ZSwgZW4gdG9kb3MgbG9zIGNhc29zLCB0ZW5lbW9zIHVuIGVycm9yIGFzb2NpYWRvIHBhcmEgbG9zIGNhbHVjbG9zIGRlIGxvcyBjb2VmaWNpZW50ZXMgZGUgcmVncmVzacOzbi4KCkFuYWxpemFtb3MgbG9zIGVycm9yZXMgZXN0w6FuZGFyICRTRSQgcGFyYSBlbCBwcmltZXIgbW9kZWxvLgoKKiBFcnJvciBlc3RhbmRhciAkXGJldGFfMCQ6IGByIHN1bW1hcnkoc2ltX3JlZ3JfMSkkY29lZmZpY2llbnRzWzNdYAoqIEVycm9yIGVzdGFuZGFyICRcYmV0YV8xJDogYHIgc3VtbWFyeShzaW1fcmVncl8xKSRjb2VmZmljaWVudHNbNF1gCgpQb2RlbW9zIGVzdGFibGVjZXIgaW50ZXJ2YWxvcyBkZSBjb25naWFuemEgYWwgJDk1XCUkIHkgJDk5XCUkIHJlc3BlY3RpdmFtZW50ZSBwYXJhIGxvcyBjb2VmaWNpZW50ZXMgZGUgcmVncmVzacOzbgoKJCQgSUNfezk1XCV9ID0gWyBcaGF0e1xiZXRhX2l9IC0gMlNFKFxoYXR7XGJldGFfaX0pLCBcaGF0e1xiZXRhX2l9ICsgMlNFKFxoYXR7XGJldGFfaX0pIF0gJCQKJCQgSUNfezk5XCV9ID0gWyBcaGF0e1xiZXRhX2l9IC0gM1NFKFxoYXR7XGJldGFfaX0pLCBcaGF0e1xiZXRhX2l9ICsgM1NFKFxoYXR7XGJldGFfaX0pIF0gJCQKCnwgUGFyYW1ldHJvIHwgSW50ZXJ2YWxvIHwgCnwtLS0tLS0tLS0tLXwtLS0tLS0tLS0tLXwKfCRcYmV0YV8wJCB8IFtgciBzaW1fcmVncl8xJGNvZWZmaWNpZW50c1sxXSAtIDIqc3VtbWFyeShzaW1fcmVncl8xKSRjb2VmZmljaWVudHNbM11gLGByIHNpbV9yZWdyXzEkY29lZmZpY2llbnRzWzFdKzIqc3VtbWFyeShzaW1fcmVncl8xKSRjb2VmZmljaWVudHNbM11gXXwKfCRcYmV0YV8xJHwgW2ByIHNpbV9yZWdyXzEkY29lZmZpY2llbnRzWzJdLTIqc3VtbWFyeShzaW1fcmVncl8xKSRjb2VmZmljaWVudHNbNF1gLGByIHNpbV9yZWdyXzEkY29lZmZpY2llbnRzWzJdKzIqc3VtbWFyeShzaW1fcmVncl8xKSRjb2VmZmljaWVudHNbNF1gXXwKCiMjIyBSZWxhY2lvbiBlbnRyZSBYIGUgWSBlbiBNQ08KU2UgcmVhbGl6YSB1biB0ZXN0IGRlIGhpcMOzdGVzaXMgbnVsYSBwYXJhIGVsIGNvZWZpY2llbnRlIGRlIHJlZ3Jlc2nDs24gJFxiZXRhXzEkIHF1ZSByZWxhY2lvbmEgbG9zIGRvcyBwYXLDoW1ldHJvcyBkZWwgbW9kZWxvLiBFbiBlc3RlIGNhc28gcHVudHVhbCBzZXLDrWFuICRzYWxlcyQgeSAkVFYkLiAkSF8wJCBlcyBsYSBoaXDDs3Rlc2lzIG51bGEgeSAkSF9hJCBlcyBsYSBoaXDDs3Rlc2lzIGFsdGVybmF0aXZhLgoKKiAkSF8wOiQgTm8gaGF5IHJlbGFjaW9uIGVudHJlICRzYWxlcyQgeSAkVFYkCiogJEhfYTokIEhheSByZWxhY2lvbiBlbnRyZSAkc2FsZXMkIHkgJFRWJAoKU2kgJFxiZXRhXzEgPSAwJCBlbnRvbmNlcyBubyBleGlzdGUgcmVsYWNpb24gZW50cmUgYW1ib3MgcGFyw6FtZXRyb3MuCgoqICRIXzA6IFxiZXRhXzEgPSAwJAoqICRIX2E6IFxiZXRhXzEgXG5lcSAwJAoKU2kgJFxiZXRhXzEkIGVzIGVzdGFkw61zdGljYW1lbnRlIHNpZ25pZmljYXRpdm8gZW50b25jZXMgZGVjaW1vcyBxdWUgZXhwbGljYSBhIGxhIHZhcmlhbnphIGRlIGxhIHZpYXJhYmxlIGRlcGVuZGllbnRlICRzYWxlcyQuCgp8IFBhcmFtZXRybyB8IFAtdmFsb3IgfCAKfC0tLS0tLS0tLS0tfC0tLS0tLS0tLS0tfAp8JFxiZXRhXzAkfCBgciBzdW1tYXJ5KHNpbV9yZWdyXzEpJGNvZWZmaWNpZW50c1s3XWB8CnwkXGJldGFfMSR8IGByIHN1bW1hcnkoc2ltX3JlZ3JfMSkkY29lZmZpY2llbnRzWzhdYHwKCiMjIFJlZ3Jlc2lvbiBMaW5lYWwgTXVsdGlwbGUKClBvZGVtb3MgYW5hbGl6YXIgZ3LDoWZpY2FtZW50ZSBsYSByZWxhY2nDs24gZW50cmUgbGFzIGRpc3RpbnRhcyB2YXJpYWJsZXMgZGUgbGEgYmFzZS4KCmBgYHtyfQpwbG90KGRhdGEpCmBgYAoKYGBge3J9CnJsbSA8LSBsbShzYWxlcyB+IHJhZGlvICsgVFYgKyBuZXdzcGFwZXIsIGRhdGE9ZGF0YSkKc3VtbWFyeShybG0pCmBgYAoKJFx0ZXh0YmZ7Q29uY2x1c2lvbmVzfSQ6IExvcyBjb2VmaWNpZW50ZXMgdGllbmVuIGxhIG1pc21hIHVuaWRhZCBxdWUgbGEgdmFyaWFibGUgaW5kZXBlbmRpZW50ZSwgeWEgc2VhIGVzdGEgJFwlJCwgJFwkJCwgdSBvdHJhLiBTZWd1biBsYSByZWdyZXNpw7NuIHJlYWxpemFkYSwgdW5hIGludmVyc2nDs24gZGUgMSB1bmlkYWQgZW4gcmFkaW8sIGVzIGxvIHF1ZSBnZW5lcmFyw6EgbcOhcyBhdW1lbnRvIGVuICRzYWxlcyQuIFRhbWJpZW4gcG9kZW1vcyBxdWUgJG5ld3NwYXBlciQgbm8gZXMgc2lnbmlmaWNhdGl2byBlc3RhZMOtc3RpY2FtZW50ZSBlbiB0w6lybWlub3MgZGUgbGEgdmFyaWFjacOzbiBkZSAkc2FsZXMkLgoKUmVhbGl6YW1vcyBsYSByZWdyZXNpb24gbGluZWFsIHF1aXRhbmRvIGEgbGEgdmFyaWFibGUgJG5ld3NwYXBlciQgcG9yIHNlciBlc3RhZMOtc3RpY2FtZW50ZSBwb2NvIHNpZ25pZmljYXRpdmEuIFNlIG9ic2VydmEgcXVlIGxvcyBjb2VmaWNpZW50ZXMgbm8gc3VmcmVuIHZhcmlhY2lvbmVzIGltcG9ydGFudGVzLgoKYGBge3J9CnJsbTIgPC0gbG0oc2FsZXMgfiByYWRpbyArIFRWICwgZGF0YT1kYXRhKQpzdW1tYXJ5KHJsbTIpCmBgYAoKIyMgUmVncmVzaW9uIExpbmVhbCBjb24gdmFyaWFibGVzIGRpY290w7NtaWNhcwoKYGBge3J9CmNyZWRpdCA8LSByZWFkLmNzdignQ3JlZGl0LmNzdicpCgpuYW1lcyhjcmVkaXQpIDwtIHRvbG93ZXIobmFtZXMoY3JlZGl0KSkKY3JlZGl0IDwtIGNyZWRpdFssLTFdCmhlYWQoY3JlZGl0KQpgYGAKCgpTZSBtb2RlbGEgZWwgYmFsYW5jZSByZXNwZWN0byBkZSBsYSB2YXJpYWJsZSBkaWNvdG9taWNhICRcdGV4dGJme2dlbmRlcn0kLiBMYSBjYXRlZ29yw61hIGJhc2Ugc2Vyw6EgJFx0ZXh0YmZ7ZmVtYWxlfSQuCmBgYHtyfQpmaXQ0IDwtIGxtKGJhbGFuY2UgfiBnZW5kZXIsIGRhdGE9Y3JlZGl0KQpzdW1tYXJ5KGZpdDQpJGNvZWZmaWNpZW50cwpgYGAKClBvZGVtb3MgaW50ZXJwcmV0YXIgYSAkXGJldGFfMCQgY29tbyBlbCBiYWxhbmNlIGRlbCBzdWJtb2RlbG8gcXVlIHJlcHJlc2VudGEgYSBsYXMgbXVqZXJlcywgeSAkXGJldGFfMSQgY29tbyBsYSB2YXJpYW56YSBxdWUgdGllbmUgZWwgYmFsYW5jZSBkZSB1bmEgcGVyc29uYSBwb3Igc2VyIGhvbWJyZS4KCiQkCmdlbmRlcl9pID0gXGxlZnRce1xiZWdpbnttYXRyaXh9CjAgJiBGZW1hbGVcXAoxICYgTWFsZQpcZW5ke21hdHJpeH1ccmlnaHQuJCQKJCQKYmFsYW5jZSA9IFxiZXRhXzAgKyBcYmV0YV8xIFxjZG90IGdlbmRlciA9IFxsZWZ0XHtcYmVnaW57bWF0cml4fQpcYmV0YV8wICsgXGJldGFfMSAmIG1hbGVcXApcYmV0YV8xICYgZmVtYWxlClxlbmR7bWF0cml4fVxyaWdodC4KJCQKClkgc2kgcXVlcmVtb3MgY2FtYmlhciBsYSBjYXRlZ29yw61hIGJhc2UgbyBsb3MgdmFsb3JlcyBkZSBsYXMgdmFyaWFibGVzIGRpY290w7NtaWNhcz8KCiQkCmdlbmRlcl9pID0gXGxlZnRce1xiZWdpbnttYXRyaXh9Ci0xICYgZmVtYWxlXFwKMSAmIG1hbGUKXGVuZHttYXRyaXh9XHJpZ2h0LgokJAoKRXN0byBkZWphIHVuYSBudWV2YSBleHByZXNpw7NuIGRlbCBtb2RlbG86CiQkCmJhbGFuY2UgPSBcYmV0YV8wICsgXGJldGFfMSBcY2RvdCBnZW5kZXIgPSBcbGVmdFx7XGJlZ2lue21hdHJpeH0KXGJldGFfMCArIFxiZXRhXzEgJiBtYWxlXFwKXGJldGFfMCAtIFxiZXRhXzEgJiBmZW1hbGUKXGVuZHttYXRyaXh9XHJpZ2h0LgokJAoKYGBge3J9CmNyZWRpdCRjb2RlZF9nZW5kZXIgPC0gaWZlbHNlKGNyZWRpdCRnZW5kZXIgPT0gIkZlbWFsZSIsLTEsMSkKZml0NSA8LSBsbShiYWxhbmNlfmNvZGVkX2dlbmRlciwgZGF0YT1jcmVkaXQpCnN1bW1hcnkoZml0NSkkY29lZmZpY2llbnRzCmBgYAoKRW4gZXN0ZSBjYXNvLCBsb3MgdmFsb3JlcyBkZSBsb3MgY29lZmljaWVudGVzIG5vIGNhbWJpYW4sIHNvbG8gc3UgaW50ZXJwcmV0YWNpw7NuLiBZYSBxdWUgJFxiZXRhXzEkIGVuIGVzdGUgY2FzbyBzZXLDoSB2YXJpYW56YSBwb3NpdGljYSBwYXJhIGfDqW5lcm8gbXVqZXIgeSB2YXJpYW56YSBuZWdhdGl2YSBwYXJhIGfDqW5lcm8gaG9tYnJlLgoKIyMgVmFyaWFibGVzIE11bHRpbm9taWFsZXMKCkVzdGUgdGlwbyBkZSB2YXJpYWJsZXMgZGlzY3JldGFzLCBxdWUgcHVlZGVuIHRvbWFyIGZpbml0b3MgdmFsb3Jlcywgc2UgZGljb3RvbWl6YW4uCgpFamVtcGxvOiBUZW5lbW9zIHVuYSB2YXJpYWJsZSBsbGFtYWRhICRcdGV4dGl0e2V0aG5pdGljeX0kIGNvbiAzIGNhdGVnb3LDrWFzLgoKYGBge3J9CiNPYnRlbmVtb3MgdmFsb3JlcyBxdWUgcHVlZGUgdG9tYXIgbGEgdmFyaWFibGUKdW5pcXVlKGNyZWRpdCRldGhuaWNpdHkpCmBgYAoKJCQKXGJldGFfMSA9IGFzaWFuX2kgPSBcbGVmdFx7XGJlZ2lue21hdHJpeH0KMCAmIG5vXFwKMSAmIHllcwpcZW5ke21hdHJpeH1ccmlnaHQuJCQKCiQkClxiZXRhXzIgPSBDYXVjYXNpYW5faSA9IFxsZWZ0XHtcYmVnaW57bWF0cml4fQowICYgbm9cXAoxICYgeWVzClxlbmR7bWF0cml4fVxyaWdodC4KJCQKCmBgYHtyfQpjcmVkaXQkY2F1Y2FzaWFuIDwtIGlmZWxzZShjcmVkaXQkZXRobmljaXR5ID09ICJDYXVjYXNpYW4iLDEsMCkKY3JlZGl0JGFzaWFuIDwtIGlmZWxzZShjcmVkaXQkZXRobmljaXR5ID09ICJBc2lhbiIsMSwwKQpgYGAKCiQkCmJhbGFuY2UgPSBcYmV0YV8wICsgXGJldGFfMSBcY2RvdCBhc2lhbiArIFxiZXRhXzIgXGNkb3QgY2F1Y2FzaWFuID0gXGxlZnRce1xiZWdpbnttYXRyaXh9ClxiZXRhXzAgKyBcYmV0YV8xICYgYXNpYW5cXApcYmV0YV8wICsgXGJldGFfMiAmIGNhdWNhc2lhblxcClxiZXRhXzAgJiBhZnJpY2FuXF9hbWVyaWNhbgpcZW5ke21hdHJpeH1ccmlnaHQuCiQkCgpgYGB7cn0KZml0NiA8LSBsbShiYWxhbmNlIH4gYXNpYW4gKyBjYXVjYXNpYW4sIGRhdGE9Y3JlZGl0KQpzdW1tYXJ5KGZpdDYpCmBgYAoK