Ejercicios Seccion 4.7

10) Esta pregunta debe responderse utilizando el conjunto de datos Semanales, que Es parte del paquete ISLR. Estos datos son similares en naturaleza a los datos de Smarket del laboratorio de este capítulo, excepto que contiene 1,089 retornos semanales de 21 años, desde principios de 1990 hasta finales de 2010.

a) Producir algunos resúmenes numéricos y gráficos de los datos semanales. ¿Parece que hay algún patrón?

library(ISLR)
library(MASS)
str(Weekly)
'data.frame':   1089 obs. of  9 variables:
 $ Year     : num  1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
 $ Lag1     : num  0.816 -0.27 -2.576 3.514 0.712 ...
 $ Lag2     : num  1.572 0.816 -0.27 -2.576 3.514 ...
 $ Lag3     : num  -3.936 1.572 0.816 -0.27 -2.576 ...
 $ Lag4     : num  -0.229 -3.936 1.572 0.816 -0.27 ...
 $ Lag5     : num  -3.484 -0.229 -3.936 1.572 0.816 ...
 $ Volume   : num  0.155 0.149 0.16 0.162 0.154 ...
 $ Today    : num  -0.27 -2.576 3.514 0.712 1.178 ...
 $ Direction: Factor w/ 2 levels "Down","Up": 1 1 2 2 2 1 2 2 2 1 ...
summary(Weekly)
      Year           Lag1               Lag2               Lag3               Lag4         
 Min.   :1990   Min.   :-18.1950   Min.   :-18.1950   Min.   :-18.1950   Min.   :-18.1950  
 1st Qu.:1995   1st Qu.: -1.1540   1st Qu.: -1.1540   1st Qu.: -1.1580   1st Qu.: -1.1580  
 Median :2000   Median :  0.2410   Median :  0.2410   Median :  0.2410   Median :  0.2380  
 Mean   :2000   Mean   :  0.1506   Mean   :  0.1511   Mean   :  0.1472   Mean   :  0.1458  
 3rd Qu.:2005   3rd Qu.:  1.4050   3rd Qu.:  1.4090   3rd Qu.:  1.4090   3rd Qu.:  1.4090  
 Max.   :2010   Max.   : 12.0260   Max.   : 12.0260   Max.   : 12.0260   Max.   : 12.0260  
      Lag5              Volume            Today          Direction 
 Min.   :-18.1950   Min.   :0.08747   Min.   :-18.1950   Down:484  
 1st Qu.: -1.1660   1st Qu.:0.33202   1st Qu.: -1.1540   Up  :605  
 Median :  0.2340   Median :1.00268   Median :  0.2410             
 Mean   :  0.1399   Mean   :1.57462   Mean   :  0.1499             
 3rd Qu.:  1.4050   3rd Qu.:2.05373   3rd Qu.:  1.4050             
 Max.   : 12.0260   Max.   :9.32821   Max.   : 12.0260             
plot(Weekly)
cor(Weekly[,-9])
              Year         Lag1        Lag2        Lag3         Lag4         Lag5
Year    1.00000000 -0.032289274 -0.03339001 -0.03000649 -0.031127923 -0.030519101
Lag1   -0.03228927  1.000000000 -0.07485305  0.05863568 -0.071273876 -0.008183096
Lag2   -0.03339001 -0.074853051  1.00000000 -0.07572091  0.058381535 -0.072499482
Lag3   -0.03000649  0.058635682 -0.07572091  1.00000000 -0.075395865  0.060657175
Lag4   -0.03112792 -0.071273876  0.05838153 -0.07539587  1.000000000 -0.075675027
Lag5   -0.03051910 -0.008183096 -0.07249948  0.06065717 -0.075675027  1.000000000
Volume  0.84194162 -0.064951313 -0.08551314 -0.06928771 -0.061074617 -0.058517414
Today  -0.03245989 -0.075031842  0.05916672 -0.07124364 -0.007825873  0.011012698
            Volume        Today
Year    0.84194162 -0.032459894
Lag1   -0.06495131 -0.075031842
Lag2   -0.08551314  0.059166717
Lag3   -0.06928771 -0.071243639
Lag4   -0.06107462 -0.007825873
Lag5   -0.05851741  0.011012698
Volume  1.00000000 -0.033077783
Today  -0.03307778  1.000000000
attach(Weekly)
The following objects are masked from Weekly (pos = 4):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

The following objects are masked from Weekly (pos = 5):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year
plot(Volume)
plot(Year)

En base a lo anterior se puede decir que el Volumen y Año tienen alguna relacion.

b) Utilice el conjunto de datos completo para realizar una regresión logística con Direction como respuesta y las cinco variables de retraso más el Volume como predictores. Utilice la función summary para imprimir los resultados. ¿Alguno de los predictores parece ser estadísticamente significativo? ¿De ser asi, cuales?

attach(Weekly)
The following objects are masked from Weekly (pos = 3):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

The following objects are masked from Weekly (pos = 5):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year

The following objects are masked from Weekly (pos = 6):

    Direction, Lag1, Lag2, Lag3, Lag4, Lag5, Today, Volume, Year
fit.glm <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Weekly, 
    family = binomial)
summary(fit.glm)

Call:
glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
    Volume, family = binomial, data = Weekly)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6949  -1.2565   0.9913   1.0849   1.4579  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  0.26686    0.08593   3.106   0.0019 **
Lag1        -0.04127    0.02641  -1.563   0.1181   
Lag2         0.05844    0.02686   2.175   0.0296 * 
Lag3        -0.01606    0.02666  -0.602   0.5469   
Lag4        -0.02779    0.02646  -1.050   0.2937   
Lag5        -0.01447    0.02638  -0.549   0.5833   
Volume      -0.02274    0.03690  -0.616   0.5377   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1496.2  on 1088  degrees of freedom
Residual deviance: 1486.4  on 1082  degrees of freedom
AIC: 1500.4

Number of Fisher Scoring iterations: 4

De los predictores utilizado el unico estadisticamente significativo es Lag2, el cual tiene un valor de P menor a 0.05

c) Calcular la matriz de confusión y la fracción global de predicciones correctas. Explique qué le está diciendo la matriz de confusión sobre los tipos de errores cometidos por la regresión logística.

probs <- predict(fit.glm, type = "response")
pred.glm <- rep("Down", length(probs))
pred.glm[probs > 0.5] <- "Up"
table(pred.glm, Direction)
        Direction
pred.glm Down  Up
    Down   54  48
    Up    430 557
mean(pred.glm == Weekly$Direction )
[1] 0.5610652
c((557/(557+48))*100,(54/(430+54))*100 ) 
[1] 92.06612 11.15702

Cuando el mercado sube la regresion logistica tiene un porcentaje de certeza del 92.1% y cuando el mercado baja el modelo funciona un 11.15%

d) Ahora ajuste el modelo de regresión logística utilizando un período de datos de entrenamiento de 1990 a 2008, con Lag2 como el único predictor. Calcule la matriz de confusión y la fracción general de las predicciones correctas para los datos retenidos (es decir, los datos de 2009 y 2010).

train <- Weekly[, "Year"] <= 2008
fit.glm2 <- glm(Direction ~ Lag2, data = Weekly, family = binomial, subset = train)
summary(fit.glm2)

Call:
glm(formula = Direction ~ Lag2, family = binomial, data = Weekly, 
    subset = train)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.536  -1.264   1.021   1.091   1.368  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  0.20326    0.06428   3.162  0.00157 **
Lag2         0.05810    0.02870   2.024  0.04298 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1354.7  on 984  degrees of freedom
Residual deviance: 1350.5  on 983  degrees of freedom
AIC: 1354.5

Number of Fisher Scoring iterations: 4
train <- (Year < 2009)
Weekly.2009_2010 <- Weekly[!train, ]
Direction.2009_2010 <- Direction[!train]
glm.probs <- predict(fit.glm2, Weekly.2009_2010,type = "response")
pred.glm2 <- rep("Down", length(glm.probs))
pred.glm2[glm.probs > 0.5] <- "Up"
table(pred.glm2,Direction.2009_2010)
         Direction.2009_2010
pred.glm2 Down Up
     Down    9  5
     Up     34 56
c(((9+56)/104)*100, (56/(56+5))*100)
[1] 62.50000 91.80328

Se tiene un 62.5% de predecciones correctas en los datos de prueba y cuando el mercado sube el modelo tiene un 91.8% de efectividad.

e) Repita (d) utilizando LDA.

fit.lda <- lda(Direction ~ Lag2, data = Weekly, subset = train)
fit.lda
Call:
lda(Direction ~ Lag2, data = Weekly, subset = train)

Prior probabilities of groups:
     Down        Up 
0.4477157 0.5522843 

Group means:
            Lag2
Down -0.03568254
Up    0.26036581

Coefficients of linear discriminants:
           LD1
Lag2 0.4414162
pred.lda <- predict(fit.lda, Weekly.2009_2010)
table(pred.lda$class, Direction.2009_2010)
      Direction.2009_2010
       Down Up
  Down    9  5
  Up     34 56

f) Repita (d) utilizando QDA

fit.qda <- qda(Direction ~ Lag2, data = Weekly, subset = train)
fit.qda
Call:
qda(Direction ~ Lag2, data = Weekly, subset = train)

Prior probabilities of groups:
     Down        Up 
0.4477157 0.5522843 

Group means:
            Lag2
Down -0.03568254
Up    0.26036581
pred.qda <- predict(fit.qda, Weekly.2009_2010)
table(pred.qda$class, Direction.2009_2010)
      Direction.2009_2010
       Down Up
  Down    0  0
  Up     43 61

g) Repita (d) utilizando KNN con K=1

library(class)
train.X <- as.matrix(Lag2[train])
test.X <- as.matrix(Lag2[!train])
train.Direction <- Direction[train]
set.seed(1)
pred.knn <- knn(train.X, test.X, train.Direction, k = 1)
table(pred.knn, Direction.2009_2010)
        Direction.2009_2010
pred.knn Down Up
    Down   21 30
    Up     22 31

h) ¿Cuál de estos métodos parece proporcionar los mejores resultados en estos datos?

Si nos basamos en la comparacion de las tasas de error de prueba, observamos que la regresión logística y LDA tienen las tasas de error mínimas, seguidas de QDA y KNN.

i) Experimente con diferentes combinaciones de predictores, incluidas posibles transformaciones e interacciones, para cada uno de los métodos. Informe las variables, el método y la matriz de confusión asociada que parece proporcionar los mejores resultados en los datos retenidos. Tenga en cuenta que también debe experimentar con valores para K en el clasificador KNN.

# Regresion logistica con Lag2:Lag1
fit.glm3 <- glm(Direction ~ Lag2:Lag1, data = Weekly, family = binomial, subset = train)
probs3 <- predict(fit.glm3, Weekly.2009_2010, type = "response")
pred.glm3 <- rep("Down", length(probs3))
pred.glm3[probs3 > 0.5] = "Up"
table(pred.glm3, Direction.2009_2010)
         Direction.2009_2010
pred.glm3 Down Up
     Down    1  1
     Up     42 60
mean(pred.glm3 == Direction.2009_2010)
[1] 0.5865385
# LDA  Lag2 con interacion  Lag1
fit.lda2 <- lda(Direction ~ Lag2:Lag1, data = Weekly, subset = train)
pred.lda2 <- predict(fit.lda2, Weekly.2009_2010)
mean(pred.lda2$class == Direction.2009_2010)
[1] 0.5769231
# QDA con sqrt(abs(Lag2))
fit.qda2 <- qda(Direction ~ Lag2 + sqrt(abs(Lag2)), data = Weekly, subset = train)
pred.qda2 <- predict(fit.qda2, Weekly.2009_2010)
table(pred.qda2$class, Direction.2009_2010)
      Direction.2009_2010
       Down Up
  Down   12 13
  Up     31 48
mean(pred.qda2$class == Direction.2009_2010)
[1] 0.5769231
# KNN k =100
pred.knn2 <- knn(train.X, test.X, train.Direction, k = 100)
table(pred.knn2, Direction.2009_2010)
         Direction.2009_2010
pred.knn2 Down Up
     Down    9 12
     Up     34 49
mean(pred.knn2 == Direction.2009_2010)
[1] 0.5576923

11) En este problema, desarrollará un modelo para predecir si un automóvil determinado obtiene un kilometraje de gasolina alto o bajo en función del conjunto de datos Automático.

a) Cree una variable binaria, mpg01, que contenga un 1 si mpg contiene un valor por encima de su mediana, y un 0 si mpg contiene un valor por debajo de su mediana. Puedes calcular la mediana usando la función median (). Tenga en cuenta que puede resultarle útil utilizar la función data.frame () para crear un único conjunto de datos que contenga mpg01 y las otras variables automáticas.

table(mpg01)
mpg01
  0   1 
196 196 

b) Explore los datos gráficamente para investigar la asociación entre mpg01 y las otras características. ¿Cuál de las otras características parece más útil para predecir mpg01? Los diagramas de dispersión y los diagramas de caja pueden ser herramientas útiles para responder esta pregunta. Describe tus hallazgos

cor(Auto[, -9])
                    mpg  cylinders displacement horsepower     weight acceleration
mpg           1.0000000 -0.7776175   -0.8051269 -0.7784268 -0.8322442    0.4233285
cylinders    -0.7776175  1.0000000    0.9508233  0.8429834  0.8975273   -0.5046834
displacement -0.8051269  0.9508233    1.0000000  0.8972570  0.9329944   -0.5438005
horsepower   -0.7784268  0.8429834    0.8972570  1.0000000  0.8645377   -0.6891955
weight       -0.8322442  0.8975273    0.9329944  0.8645377  1.0000000   -0.4168392
acceleration  0.4233285 -0.5046834   -0.5438005 -0.6891955 -0.4168392    1.0000000
year          0.5805410 -0.3456474   -0.3698552 -0.4163615 -0.3091199    0.2903161
origin        0.5652088 -0.5689316   -0.6145351 -0.4551715 -0.5850054    0.2127458
mpg01         0.8369392 -0.7591939   -0.7534766 -0.6670526 -0.7577566    0.3468215
mpg01.1       0.8369392 -0.7591939   -0.7534766 -0.6670526 -0.7577566    0.3468215
mpg01.2       0.8369392 -0.7591939   -0.7534766 -0.6670526 -0.7577566    0.3468215
mpg01.3       0.8369392 -0.7591939   -0.7534766 -0.6670526 -0.7577566    0.3468215
                   year     origin      mpg01    mpg01.1    mpg01.2    mpg01.3
mpg           0.5805410  0.5652088  0.8369392  0.8369392  0.8369392  0.8369392
cylinders    -0.3456474 -0.5689316 -0.7591939 -0.7591939 -0.7591939 -0.7591939
displacement -0.3698552 -0.6145351 -0.7534766 -0.7534766 -0.7534766 -0.7534766
horsepower   -0.4163615 -0.4551715 -0.6670526 -0.6670526 -0.6670526 -0.6670526
weight       -0.3091199 -0.5850054 -0.7577566 -0.7577566 -0.7577566 -0.7577566
acceleration  0.2903161  0.2127458  0.3468215  0.3468215  0.3468215  0.3468215
year          1.0000000  0.1815277  0.4299042  0.4299042  0.4299042  0.4299042
origin        0.1815277  1.0000000  0.5136984  0.5136984  0.5136984  0.5136984
mpg01         0.4299042  0.5136984  1.0000000  1.0000000  1.0000000  1.0000000
mpg01.1       0.4299042  0.5136984  1.0000000  1.0000000  1.0000000  1.0000000
mpg01.2       0.4299042  0.5136984  1.0000000  1.0000000  1.0000000  1.0000000
mpg01.3       0.4299042  0.5136984  1.0000000  1.0000000  1.0000000  1.0000000
par(mfrow = c(2,2))
boxplot(acceleration~mpg01,Data,main="Accel~mpg01")
boxplot(weight~mpg01,Data,main="Weight~mpg01")
boxplot(horsepower~mpg01,Data,main="Horsep~mpg01")
boxplot(displacement~mpg01,Data,main="Disp~mpg01")

boxplot(cylinders~mpg01,Data,main="Disp~mpg01")
boxplot(year~mpg01,Data,main="Disp~mpg01")

Se puede observar una asociacion entre mpg01 y cylinders, weight, displacement y horsepower

c) Dividir los datos en un conjunto de entrenamiento y un conjunto de prueba

train <- (year %% 2 == 0)
Auto.train <- Auto[train, ]
Auto.test <- Auto[!train, ]
mpg01.test <- mpg01[!train]

d) Realice LDA en los datos de entrenamiento para predecir “mpg01” utilizando las variables que parecían estar más asociadas con “mpg01” en (b). ¿Cuál es el error de prueba del modelo obtenido?

fit.lda <- lda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, 
               subset = train)
fit.lda
Call:
lda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, 
    subset = train)

Prior probabilities of groups:
        0         1 
0.4571429 0.5428571 

Group means:
  cylinders   weight displacement horsepower
0  6.812500 3604.823     271.7396  133.14583
1  4.070175 2314.763     111.6623   77.92105

Coefficients of linear discriminants:
                       LD1
cylinders    -0.6741402638
weight       -0.0011465750
displacement  0.0004481325
horsepower    0.0059035377
pred.lda <- predict(fit.lda, Auto.test)
table(pred.lda$class, mpg01.test)
   mpg01.test
     0  1
  0 86  9
  1 14 73
(mean(pred.lda$class != mpg01.test)*100)
[1] 12.63736

la tasa de error de prueba es del 12,6%

e) Realice QDA en los datos de entrenamiento para predecir “mpg01” utilizando las variables que parecían estar más asociadas con “mpg01” en (b). ¿Cuál es el error de prueba del modelo obtenido?

fit.qda <- qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto,
               subset = train)
fit.qda
Call:
qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, 
    subset = train)

Prior probabilities of groups:
        0         1 
0.4571429 0.5428571 

Group means:
  cylinders   weight displacement horsepower
0  6.812500 3604.823     271.7396  133.14583
1  4.070175 2314.763     111.6623   77.92105
pred.qda <- predict(fit.qda, Auto.test)
table(pred.qda$class, mpg01.test)
   mpg01.test
     0  1
  0 89 13
  1 11 69
(mean(pred.qda$class != mpg01.test)*100)
[1] 13.18681

la tasa de error de prueba es del 13,18%

f) Realice una regresión logística en los datos de entrenamiento para predecir “mpg01” utilizando las variables que parecían estar más asociadas con “mpg01” en (b). ¿Cuál es el error de prueba del modelo obtenido?

fit.glm <- glm(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, family = binomial, subset = train)
summary(fit.glm)

Call:
glm(formula = mpg01 ~ cylinders + weight + displacement + horsepower, 
    family = binomial, data = Auto, subset = train)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.48027  -0.03413   0.10583   0.29634   2.57584  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  17.658730   3.409012   5.180 2.22e-07 ***
cylinders    -1.028032   0.653607  -1.573   0.1158    
weight       -0.002922   0.001137  -2.569   0.0102 *  
displacement  0.002462   0.015030   0.164   0.8699    
horsepower   -0.050611   0.025209  -2.008   0.0447 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 289.58  on 209  degrees of freedom
Residual deviance:  83.24  on 205  degrees of freedom
AIC: 93.24

Number of Fisher Scoring iterations: 7
probs <- predict(fit.glm, Auto.test, type = "response")
pred.glm <- rep(0, length(probs))
pred.glm[probs > 0.5] <- 1
table(pred.glm, mpg01.test)
        mpg01.test
pred.glm  0  1
       0 89 11
       1 11 71
(mean(pred.glm != mpg01.test)*100)
[1] 12.08791

La tasa de error es del 12,1%

g) Ejecute KNN en los datos de entrenamiento, con varios valores de K, para predecir “mpg01” utilizando las variables que parecían estar más asociadas con “mpg01” en (b). ¿Qué errores de prueba obtienes? ¿Qué valor de K parece funcionar mejor en este conjunto de datos?

train.X <- cbind(cylinders, weight, displacement, horsepower)[train, ]
test.X <- cbind(cylinders, weight, displacement, horsepower)[!train, ]
train.mpg01 <- mpg01[train]
set.seed(1)
pred.knn <- knn(train.X, test.X, train.mpg01, k = 1)
table(pred.knn, mpg01.test)
        mpg01.test
pred.knn  0  1
       0 83 11
       1 17 71
mean(pred.knn != mpg01.test)*100
[1] 15.38462
pred.knn <- knn(train.X, test.X, train.mpg01, k = 10)
table(pred.knn, mpg01.test)
        mpg01.test
pred.knn  0  1
       0 77  7
       1 23 75
mean(pred.knn != mpg01.test)*100
[1] 16.48352
pred.knn <- knn(train.X, test.X, train.mpg01, k = 100)
table(pred.knn, mpg01.test)
        mpg01.test
pred.knn  0  1
       0 81  7
       1 19 75
mean(pred.knn != mpg01.test)*100
[1] 14.28571

cuando k=100 se obtiene la menor tasa de error

12) Este problema implica funciones de escritura.

a) Escriba una función, Power (), que imprima el resultado de aumentar 2 a la 3ª potencia. En otras palabras, su función debe calcular 2^3 e imprimir los resultados.

Power <- function() {
    2^3
}
print(Power())
[1] 8

b) Cree una nueva función, Power2(), que le permita pasar cualquiera de los dos números, “x” y “a”, e imprima el valor de “x ^ a”.

Power2 <- function(x, a) {
    x^a
}
Power2(3, 8)
[1] 6561

c) Usando la función Power2 () que acaba de escribir, calcule 10^3, 8^17 y 131^3.

Power2(10, 3)
[1] 1000
Power2(8, 17)
[1] 2.2518e+15
Power2(131, 3)
[1] 2248091

d) Ahora cree una nueva función, Power3 (), que realmente devuelva el resultado “x ^ a” como un objeto R, en lugar de simplemente imprimirlo en la pantalla. Es decir, si almacena el valor “x ^ a” en un objeto llamado “resultado” dentro de su función, entonces simplemente puede devolver () este resultado

Power3 <- function(x , a) {
    result <- x^a
    return(result)
}

e) Ahora, utilizando la función Power3(), cree una gráfica de f (x) = x^3. El eje x debe mostrar un rango de enteros de 1 a 10, y el eje y debe mostrar x^2. Etiquete los ejes apropiadamente y use un título apropiado para la figura. Considere mostrar el eje x, el eje y, o ambos en la escala logarítmica.

x <- 1:10
plot(x, Power3(x, 2), log = "xy", xlab = "Log of x", ylab = "Log of x^2", main = "Log of x^2 vs Log of x")

f) Cree una función, PlotPower (), que le permita crear un gráfico de “x” contra “x ^ a” para una “a” fija y para un rango de valores de “x”.

PlotPower <- function(x, a) {
    plot(x, Power3(x, a))
}
PlotPower(1:10, 3)

13) Usando el conjunto de datos “Boston”, ajuste los modelos de clasificación para predecir si un suburbio dado tiene una tasa de criminalidad por encima o por debajo de la mediana. Explore los modelos de regresión logística, LDA y KNN utilizando varios subconjuntos de los predictores. Describe tus hallazgos.

library(MASS)
attach(Boston)
crim01 <- rep(0, length(crim))
crim01[crim > median(crim)] <- 1
Boston <- data.frame(Boston, crim01)
train <- 1:(length(crim) / 2)
test <- (length(crim) / 2 + 1):length(crim)
Boston.train <- Boston[train, ]
Boston.test <- Boston[test, ]
crim01.test <- crim01[test]
fit.glm <- glm(crim01 ~ . - crim01 - crim, data = Boston, family = binomial, subset = train)
glm.fit: fitted probabilities numerically 0 or 1 occurred
probs <- predict(fit.glm, Boston.test, type = "response")
pred.glm <- rep(0, length(probs))
pred.glm[probs > 0.5] <- 1
table(pred.glm, crim01.test)
        crim01.test
pred.glm   0   1
       0  68  24
       1  22 139
mean(pred.glm != crim01.test)*100
[1] 18.18182
fit.glm <- glm(crim01 ~ . - crim01 - crim - chas - nox, data = Boston, family = binomial, subset = train)
glm.fit: fitted probabilities numerically 0 or 1 occurred
probs <- predict(fit.glm, Boston.test, type = "response")
pred.glm <- rep(0, length(probs))
pred.glm[probs > 0.5] <- 1
table(pred.glm, crim01.test)
        crim01.test
pred.glm   0   1
       0  78  28
       1  12 135
mean(pred.glm != crim01.test)*100
[1] 15.81028
fit.lda <- lda(crim01 ~ . - crim01 - crim, data = Boston, subset = train)
pred.lda <- predict(fit.lda, Boston.test)
table(pred.lda$class, crim01.test)
   crim01.test
      0   1
  0  80  24
  1  10 139
mean(pred.lda$class != crim01.test)*100
[1] 13.43874
fit.lda <- lda(crim01 ~ . - crim01 - crim - chas - nox, data = Boston, subset = train)
pred.lda <- predict(fit.lda, Boston.test)
table(pred.lda$class, crim01.test)
   crim01.test
      0   1
  0  82  30
  1   8 133
mean(pred.lda$class != crim01.test)*100
[1] 15.01976
train.X <- cbind(zn, indus, chas, nox, rm, age, dis, rad, tax, ptratio, black, lstat, medv)[train, ]
test.X <- cbind(zn, indus, chas, nox, rm, age, dis, rad, tax, ptratio, black, lstat, medv)[test, ]
train.crim01 <- crim01[train]
set.seed(1)
pred.knn <- knn(train.X, test.X, train.crim01, k = 1)
table(pred.knn, crim01.test)
        crim01.test
pred.knn   0   1
       0  85 111
       1   5  52
mean(pred.knn != crim01.test)*100
[1] 45.8498
pred.knn <- knn(train.X, test.X, train.crim01, k = 10)
table(pred.knn, crim01.test)
        crim01.test
pred.knn   0   1
       0  83  23
       1   7 140
mean(pred.knn != crim01.test)*100
[1] 11.85771
pred.knn <- knn(train.X, test.X, train.crim01, k = 100)
table(pred.knn, crim01.test)
        crim01.test
pred.knn   0   1
       0  86 120
       1   4  43
mean(pred.knn != crim01.test)*100
[1] 49.01186
LS0tDQp0aXRsZTogIkhvamEgRGUgVHJhYmFqbyAjMyINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQojIyBFamVyY2ljaW9zIFNlY2Npb24gNC43DQoNCiMjIyMgMTApIEVzdGEgcHJlZ3VudGEgZGViZSByZXNwb25kZXJzZSB1dGlsaXphbmRvIGVsIGNvbmp1bnRvIGRlIGRhdG9zIFNlbWFuYWxlcywgcXVlIEVzIHBhcnRlIGRlbCBwYXF1ZXRlIElTTFIuIEVzdG9zIGRhdG9zIHNvbiBzaW1pbGFyZXMgZW4gbmF0dXJhbGV6YSBhIGxvcyBkYXRvcyBkZSBTbWFya2V0IGRlbCBsYWJvcmF0b3JpbyBkZSBlc3RlIGNhcMOtdHVsbywgZXhjZXB0byBxdWUgY29udGllbmUgMSwwODkgcmV0b3Jub3Mgc2VtYW5hbGVzIGRlIDIxIGHDsW9zLCBkZXNkZSBwcmluY2lwaW9zIGRlIDE5OTAgaGFzdGEgZmluYWxlcyBkZSAyMDEwLg0KDQojIyMjIGEpIFByb2R1Y2lyIGFsZ3Vub3MgcmVzw7ptZW5lcyBudW3DqXJpY29zIHkgZ3LDoWZpY29zIGRlIGxvcyBkYXRvcyBzZW1hbmFsZXMuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgwr9QYXJlY2UgcXVlIGhheSBhbGfDum4gcGF0csOzbj8NCg0KYGBge3J9DQpsaWJyYXJ5KElTTFIpDQpsaWJyYXJ5KE1BU1MpDQpzdHIoV2Vla2x5KQ0Kc3VtbWFyeShXZWVrbHkpDQpwbG90KFdlZWtseSkNCmNvcihXZWVrbHlbLC05XSkNCmF0dGFjaChXZWVrbHkpDQpwbG90KFZvbHVtZSkNCnBsb3QoWWVhcikNCg0KYGBgDQoNCkVuIGJhc2UgYSBsbyBhbnRlcmlvciBzZSBwdWVkZSBkZWNpciBxdWUgZWwgVm9sdW1lbiB5IEHDsW8gdGllbmVuIGFsZ3VuYSByZWxhY2lvbi4NCg0KDQojIyMjIGIpIFV0aWxpY2UgZWwgY29uanVudG8gZGUgZGF0b3MgY29tcGxldG8gcGFyYSByZWFsaXphciB1bmEgcmVncmVzacOzbiBsb2fDrXN0aWNhIGNvbiBEaXJlY3Rpb24gY29tbyByZXNwdWVzdGEgeSBsYXMgY2luY28gdmFyaWFibGVzIGRlIHJldHJhc28gbcOhcyBlbCBWb2x1bWUgY29tbyBwcmVkaWN0b3Jlcy4gVXRpbGljZSBsYSBmdW5jacOzbiBzdW1tYXJ5IHBhcmEgaW1wcmltaXIgbG9zIHJlc3VsdGFkb3MuIMK/QWxndW5vIGRlIGxvcyBwcmVkaWN0b3JlcyBwYXJlY2Ugc2VyIGVzdGFkw61zdGljYW1lbnRlIHNpZ25pZmljYXRpdm8/IMK/RGUgc2VyIGFzaSwgY3VhbGVzPw0KDQpgYGB7cn0NCmF0dGFjaChXZWVrbHkpDQpmaXQuZ2xtIDwtIGdsbShEaXJlY3Rpb24gfiBMYWcxICsgTGFnMiArIExhZzMgKyBMYWc0ICsgTGFnNSArIFZvbHVtZSwgZGF0YSA9IFdlZWtseSwgDQogICAgZmFtaWx5ID0gYmlub21pYWwpDQpzdW1tYXJ5KGZpdC5nbG0pDQoNCmBgYA0KDQpEZSBsb3MgcHJlZGljdG9yZXMgdXRpbGl6YWRvIGVsIHVuaWNvIGVzdGFkaXN0aWNhbWVudGUgc2lnbmlmaWNhdGl2byBlcyBMYWcyLCBlbCBjdWFsIHRpZW5lIHVuIHZhbG9yIGRlIFAgbWVub3IgYSAwLjA1DQoNCg0KIyMjIyBjKSBDYWxjdWxhciBsYSBtYXRyaXogZGUgY29uZnVzacOzbiB5IGxhIGZyYWNjacOzbiBnbG9iYWwgZGUgcHJlZGljY2lvbmVzIGNvcnJlY3Rhcy4gRXhwbGlxdWUgcXXDqSBsZSBlc3TDoSBkaWNpZW5kbyBsYSBtYXRyaXogZGUgY29uZnVzacOzbiBzb2JyZSBsb3MgdGlwb3MgZGUgZXJyb3JlcyBjb21ldGlkb3MgcG9yIGxhIHJlZ3Jlc2nDs24gbG9nw61zdGljYS4NCg0KYGBge3J9DQoNCnByb2JzIDwtIHByZWRpY3QoZml0LmdsbSwgdHlwZSA9ICJyZXNwb25zZSIpDQpwcmVkLmdsbSA8LSByZXAoIkRvd24iLCBsZW5ndGgocHJvYnMpKQ0KcHJlZC5nbG1bcHJvYnMgPiAwLjVdIDwtICJVcCINCnRhYmxlKHByZWQuZ2xtLCBEaXJlY3Rpb24pDQoNCmBgYA0KDQpgYGB7cn0NCm1lYW4ocHJlZC5nbG0gPT0gV2Vla2x5JERpcmVjdGlvbiApDQpgYGANCg0KYGBge3J9DQpjKCg1NTcvKDU1Nys0OCkpKjEwMCwoNTQvKDQzMCs1NCkpKjEwMCApIA0KDQpgYGANCg0KQ3VhbmRvIGVsIG1lcmNhZG8gc3ViZSBsYSByZWdyZXNpb24gbG9naXN0aWNhIHRpZW5lIHVuIHBvcmNlbnRhamUgZGUgY2VydGV6YSBkZWwgOTIuMSUgeSBjdWFuZG8gZWwgbWVyY2FkbyBiYWphIGVsIG1vZGVsbyBmdW5jaW9uYSB1biAxMS4xNSUNCg0KDQojIyMjIGQpIEFob3JhIGFqdXN0ZSBlbCBtb2RlbG8gZGUgcmVncmVzacOzbiBsb2fDrXN0aWNhIHV0aWxpemFuZG8gdW4gcGVyw61vZG8gZGUgZGF0b3MgZGUgZW50cmVuYW1pZW50byBkZSAxOTkwIGEgMjAwOCwgY29uIExhZzIgY29tbyBlbCDDum5pY28gcHJlZGljdG9yLiBDYWxjdWxlIGxhIG1hdHJpeiBkZSBjb25mdXNpw7NuIHkgbGEgZnJhY2Npw7NuIGdlbmVyYWwgZGUgbGFzIHByZWRpY2Npb25lcyBjb3JyZWN0YXMgcGFyYSBsb3MgZGF0b3MgcmV0ZW5pZG9zIChlcyBkZWNpciwgbG9zIGRhdG9zIGRlIDIwMDkgeSAyMDEwKS4NCg0KYGBge3J9DQp0cmFpbiA8LSBXZWVrbHlbLCAiWWVhciJdIDw9IDIwMDgNCmZpdC5nbG0yIDwtIGdsbShEaXJlY3Rpb24gfiBMYWcyLCBkYXRhID0gV2Vla2x5LCBmYW1pbHkgPSBiaW5vbWlhbCwgc3Vic2V0ID0gdHJhaW4pDQpzdW1tYXJ5KGZpdC5nbG0yKQ0KDQoNCg0KDQpgYGANCg0KYGBge3J9DQp0cmFpbiA8LSAoWWVhciA8IDIwMDkpDQpXZWVrbHkuMjAwOV8yMDEwIDwtIFdlZWtseVshdHJhaW4sIF0NCkRpcmVjdGlvbi4yMDA5XzIwMTAgPC0gRGlyZWN0aW9uWyF0cmFpbl0NCmdsbS5wcm9icyA8LSBwcmVkaWN0KGZpdC5nbG0yLCBXZWVrbHkuMjAwOV8yMDEwLHR5cGUgPSAicmVzcG9uc2UiKQ0KcHJlZC5nbG0yIDwtIHJlcCgiRG93biIsIGxlbmd0aChnbG0ucHJvYnMpKQ0KcHJlZC5nbG0yW2dsbS5wcm9icyA+IDAuNV0gPC0gIlVwIg0KDQp0YWJsZShwcmVkLmdsbTIsRGlyZWN0aW9uLjIwMDlfMjAxMCkNCmBgYA0KDQpgYGB7cn0NCmMoKCg5KzU2KS8xMDQpKjEwMCwgKDU2Lyg1Nis1KSkqMTAwKQ0KYGBgDQpTZSB0aWVuZSB1biA2Mi41JSBkZSBwcmVkZWNjaW9uZXMgY29ycmVjdGFzIGVuIGxvcyBkYXRvcyBkZSBwcnVlYmEgeSBjdWFuZG8gZWwgbWVyY2FkbyBzdWJlIGVsIG1vZGVsbyB0aWVuZSB1biA5MS44JSBkZSBlZmVjdGl2aWRhZC4NCg0KIyMjIyBlKSBSZXBpdGEgKGQpIHV0aWxpemFuZG8gTERBLg0KDQpgYGB7cn0NCmZpdC5sZGEgPC0gbGRhKERpcmVjdGlvbiB+IExhZzIsIGRhdGEgPSBXZWVrbHksIHN1YnNldCA9IHRyYWluKQ0KZml0LmxkYQ0KDQpwcmVkLmxkYSA8LSBwcmVkaWN0KGZpdC5sZGEsIFdlZWtseS4yMDA5XzIwMTApDQp0YWJsZShwcmVkLmxkYSRjbGFzcywgRGlyZWN0aW9uLjIwMDlfMjAxMCkNCg0KYGBgDQoNCg0KIyMjIyBmKSBSZXBpdGEgKGQpIHV0aWxpemFuZG8gUURBIA0KDQpgYGB7cn0NCmZpdC5xZGEgPC0gcWRhKERpcmVjdGlvbiB+IExhZzIsIGRhdGEgPSBXZWVrbHksIHN1YnNldCA9IHRyYWluKQ0KZml0LnFkYQ0KDQpwcmVkLnFkYSA8LSBwcmVkaWN0KGZpdC5xZGEsIFdlZWtseS4yMDA5XzIwMTApDQp0YWJsZShwcmVkLnFkYSRjbGFzcywgRGlyZWN0aW9uLjIwMDlfMjAxMCkNCg0KDQpgYGANCg0KIyMjIyBnKSBSZXBpdGEgKGQpIHV0aWxpemFuZG8gS05OIGNvbiBLPTENCg0KYGBge3J9DQpsaWJyYXJ5KGNsYXNzKQ0KdHJhaW4uWCA8LSBhcy5tYXRyaXgoTGFnMlt0cmFpbl0pDQp0ZXN0LlggPC0gYXMubWF0cml4KExhZzJbIXRyYWluXSkNCnRyYWluLkRpcmVjdGlvbiA8LSBEaXJlY3Rpb25bdHJhaW5dDQpzZXQuc2VlZCgxKQ0KcHJlZC5rbm4gPC0ga25uKHRyYWluLlgsIHRlc3QuWCwgdHJhaW4uRGlyZWN0aW9uLCBrID0gMSkNCnRhYmxlKHByZWQua25uLCBEaXJlY3Rpb24uMjAwOV8yMDEwKQ0KYGBgDQoNCg0KIyMjIyBoKSDCv0N1w6FsIGRlIGVzdG9zIG3DqXRvZG9zIHBhcmVjZSBwcm9wb3JjaW9uYXIgbG9zIG1lam9yZXMgcmVzdWx0YWRvcyBlbiBlc3RvcyBkYXRvcz8NCg0KU2kgbm9zIGJhc2Ftb3MgZW4gbGEgY29tcGFyYWNpb24gZGUgbGFzIHRhc2FzIGRlIGVycm9yIGRlIHBydWViYSwgb2JzZXJ2YW1vcyBxdWUgbGEgcmVncmVzacOzbiBsb2fDrXN0aWNhIHkgTERBIHRpZW5lbiBsYXMgdGFzYXMgZGUgZXJyb3IgbcOtbmltYXMsIHNlZ3VpZGFzIGRlIFFEQSB5IEtOTi4NCg0KDQojIyMjIGkpIEV4cGVyaW1lbnRlIGNvbiBkaWZlcmVudGVzIGNvbWJpbmFjaW9uZXMgZGUgcHJlZGljdG9yZXMsIGluY2x1aWRhcyBwb3NpYmxlcyB0cmFuc2Zvcm1hY2lvbmVzIGUgaW50ZXJhY2Npb25lcywgcGFyYSBjYWRhIHVubyBkZSBsb3MgbcOpdG9kb3MuIEluZm9ybWUgbGFzIHZhcmlhYmxlcywgZWwgbcOpdG9kbyB5IGxhIG1hdHJpeiBkZSBjb25mdXNpw7NuIGFzb2NpYWRhIHF1ZSBwYXJlY2UgcHJvcG9yY2lvbmFyIGxvcyBtZWpvcmVzIHJlc3VsdGFkb3MgZW4gbG9zIGRhdG9zIHJldGVuaWRvcy4gVGVuZ2EgZW4gY3VlbnRhIHF1ZSB0YW1iacOpbiBkZWJlIGV4cGVyaW1lbnRhciBjb24gdmFsb3JlcyBwYXJhIEsgZW4gZWwgY2xhc2lmaWNhZG9yIEtOTi4NCg0KYGBge3J9DQojIFJlZ3Jlc2lvbiBsb2dpc3RpY2EgY29uIExhZzI6TGFnMQ0KDQpmaXQuZ2xtMyA8LSBnbG0oRGlyZWN0aW9uIH4gTGFnMjpMYWcxLCBkYXRhID0gV2Vla2x5LCBmYW1pbHkgPSBiaW5vbWlhbCwgc3Vic2V0ID0gdHJhaW4pDQpwcm9iczMgPC0gcHJlZGljdChmaXQuZ2xtMywgV2Vla2x5LjIwMDlfMjAxMCwgdHlwZSA9ICJyZXNwb25zZSIpDQpwcmVkLmdsbTMgPC0gcmVwKCJEb3duIiwgbGVuZ3RoKHByb2JzMykpDQpwcmVkLmdsbTNbcHJvYnMzID4gMC41XSA9ICJVcCINCnRhYmxlKHByZWQuZ2xtMywgRGlyZWN0aW9uLjIwMDlfMjAxMCkNCm1lYW4ocHJlZC5nbG0zID09IERpcmVjdGlvbi4yMDA5XzIwMTApDQpgYGANCg0KYGBge3J9DQojIExEQSAgTGFnMiBjb24gaW50ZXJhY2lvbiAgTGFnMQ0KZml0LmxkYTIgPC0gbGRhKERpcmVjdGlvbiB+IExhZzI6TGFnMSwgZGF0YSA9IFdlZWtseSwgc3Vic2V0ID0gdHJhaW4pDQpwcmVkLmxkYTIgPC0gcHJlZGljdChmaXQubGRhMiwgV2Vla2x5LjIwMDlfMjAxMCkNCm1lYW4ocHJlZC5sZGEyJGNsYXNzID09IERpcmVjdGlvbi4yMDA5XzIwMTApDQpgYGANCg0KYGBge3J9DQojIFFEQSBjb24gc3FydChhYnMoTGFnMikpDQpmaXQucWRhMiA8LSBxZGEoRGlyZWN0aW9uIH4gTGFnMiArIHNxcnQoYWJzKExhZzIpKSwgZGF0YSA9IFdlZWtseSwgc3Vic2V0ID0gdHJhaW4pDQpwcmVkLnFkYTIgPC0gcHJlZGljdChmaXQucWRhMiwgV2Vla2x5LjIwMDlfMjAxMCkNCnRhYmxlKHByZWQucWRhMiRjbGFzcywgRGlyZWN0aW9uLjIwMDlfMjAxMCkNCm1lYW4ocHJlZC5xZGEyJGNsYXNzID09IERpcmVjdGlvbi4yMDA5XzIwMTApDQpgYGANCg0KYGBge3J9DQojIEtOTiBrID0xMDANCnByZWQua25uMiA8LSBrbm4odHJhaW4uWCwgdGVzdC5YLCB0cmFpbi5EaXJlY3Rpb24sIGsgPSAxMDApDQp0YWJsZShwcmVkLmtubjIsIERpcmVjdGlvbi4yMDA5XzIwMTApDQptZWFuKHByZWQua25uMiA9PSBEaXJlY3Rpb24uMjAwOV8yMDEwKQ0KYGBgDQoNCg0KIyMjIyAxMSkgRW4gZXN0ZSBwcm9ibGVtYSwgZGVzYXJyb2xsYXLDoSB1biBtb2RlbG8gcGFyYSBwcmVkZWNpciBzaSB1biBhdXRvbcOzdmlsIGRldGVybWluYWRvIG9idGllbmUgdW4ga2lsb21ldHJhamUgZGUgZ2Fzb2xpbmEgYWx0byBvIGJham8gZW4gZnVuY2nDs24gZGVsIGNvbmp1bnRvIGRlIGRhdG9zIEF1dG9tw6F0aWNvLg0KDQojIyMjIGEpIENyZWUgdW5hIHZhcmlhYmxlIGJpbmFyaWEsIG1wZzAxLCBxdWUgY29udGVuZ2EgdW4gMSBzaSBtcGcgY29udGllbmUgdW4gdmFsb3IgcG9yIGVuY2ltYSBkZSBzdSBtZWRpYW5hLCB5IHVuIDAgc2kgbXBnIGNvbnRpZW5lIHVuIHZhbG9yIHBvciBkZWJham8gZGUgc3UgbWVkaWFuYS4gUHVlZGVzIGNhbGN1bGFyIGxhIG1lZGlhbmEgdXNhbmRvIGxhIGZ1bmNpw7NuIG1lZGlhbiAoKS4gVGVuZ2EgZW4gY3VlbnRhIHF1ZSBwdWVkZSByZXN1bHRhcmxlIMO6dGlsIHV0aWxpemFyIGxhIGZ1bmNpw7NuIGRhdGEuZnJhbWUgKCkgcGFyYSBjcmVhciB1biDDum5pY28gY29uanVudG8gZGUgZGF0b3MgcXVlIGNvbnRlbmdhIG1wZzAxIHkgbGFzIG90cmFzIHZhcmlhYmxlcyBhdXRvbcOhdGljYXMuDQoNCmBgYHtyfQ0KYXR0YWNoKEF1dG8pDQptcGcwMSA8LSByZXAoMCwgbGVuZ3RoKG1wZykpDQptcGcwMVttcGcgPiBtZWRpYW4obXBnKV0gPC0gMQ0KQXV0byA8LSBkYXRhLmZyYW1lKEF1dG8sIG1wZzAxKQ0KQXV0bw0KdGFibGUobXBnMDEpDQpgYGANCg0KIyMjIyBiKSBFeHBsb3JlIGxvcyBkYXRvcyBncsOhZmljYW1lbnRlIHBhcmEgaW52ZXN0aWdhciBsYSBhc29jaWFjacOzbiBlbnRyZSBtcGcwMSB5IGxhcyBvdHJhcyBjYXJhY3RlcsOtc3RpY2FzLiDCv0N1w6FsIGRlIGxhcyBvdHJhcyBjYXJhY3RlcsOtc3RpY2FzIHBhcmVjZSBtw6FzIMO6dGlsIHBhcmEgcHJlZGVjaXIgbXBnMDE/IExvcyBkaWFncmFtYXMgZGUgZGlzcGVyc2nDs24geSBsb3MgZGlhZ3JhbWFzIGRlIGNhamEgcHVlZGVuIHNlciBoZXJyYW1pZW50YXMgw7p0aWxlcyBwYXJhIHJlc3BvbmRlciBlc3RhIHByZWd1bnRhLiBEZXNjcmliZSB0dXMgaGFsbGF6Z29zIA0KDQpgYGB7cn0NCmNvcihBdXRvWywgLTldKQ0KcGFyKG1mcm93ID0gYygyLDIpKQ0KYm94cGxvdChhY2NlbGVyYXRpb25+bXBnMDEsRGF0YSxtYWluPSJBY2NlbH5tcGcwMSIpDQpib3hwbG90KHdlaWdodH5tcGcwMSxEYXRhLG1haW49IldlaWdodH5tcGcwMSIpDQpib3hwbG90KGhvcnNlcG93ZXJ+bXBnMDEsRGF0YSxtYWluPSJIb3JzZXB+bXBnMDEiKQ0KYm94cGxvdChkaXNwbGFjZW1lbnR+bXBnMDEsRGF0YSxtYWluPSJEaXNwfm1wZzAxIikNCmJveHBsb3QoY3lsaW5kZXJzfm1wZzAxLERhdGEsbWFpbj0iRGlzcH5tcGcwMSIpDQpib3hwbG90KHllYXJ+bXBnMDEsRGF0YSxtYWluPSJEaXNwfm1wZzAxIikNCg0KYGBgDQoNClNlIHB1ZWRlIG9ic2VydmFyIHVuYSBhc29jaWFjaW9uIGVudHJlIG1wZzAxIHkgY3lsaW5kZXJzLCB3ZWlnaHQsIGRpc3BsYWNlbWVudCB5IGhvcnNlcG93ZXINCg0KIyMjIyBjKSBEaXZpZGlyIGxvcyBkYXRvcyBlbiB1biBjb25qdW50byBkZSBlbnRyZW5hbWllbnRvIHkgdW4gY29uanVudG8gZGUgcHJ1ZWJhDQoNCmBgYHtyfQ0KdHJhaW4gPC0gKHllYXIgJSUgMiA9PSAwKQ0KQXV0by50cmFpbiA8LSBBdXRvW3RyYWluLCBdDQpBdXRvLnRlc3QgPC0gQXV0b1shdHJhaW4sIF0NCm1wZzAxLnRlc3QgPC0gbXBnMDFbIXRyYWluXQ0KYGBgDQoNCg0KIyMjIyBkKSBSZWFsaWNlIExEQSBlbiBsb3MgZGF0b3MgZGUgZW50cmVuYW1pZW50byBwYXJhIHByZWRlY2lyICJtcGcwMSIgdXRpbGl6YW5kbyBsYXMgdmFyaWFibGVzIHF1ZSBwYXJlY8OtYW4gZXN0YXIgbcOhcyBhc29jaWFkYXMgY29uICJtcGcwMSIgZW4gKGIpLiDCv0N1w6FsIGVzIGVsIGVycm9yIGRlIHBydWViYSBkZWwgbW9kZWxvIG9idGVuaWRvPw0KDQpgYGB7cn0NCmZpdC5sZGEgPC0gbGRhKG1wZzAxIH4gY3lsaW5kZXJzICsgd2VpZ2h0ICsgZGlzcGxhY2VtZW50ICsgaG9yc2Vwb3dlciwgZGF0YSA9IEF1dG8sIA0KICAgICAgICAgICAgICAgc3Vic2V0ID0gdHJhaW4pDQpmaXQubGRhDQpwcmVkLmxkYSA8LSBwcmVkaWN0KGZpdC5sZGEsIEF1dG8udGVzdCkNCnRhYmxlKHByZWQubGRhJGNsYXNzLCBtcGcwMS50ZXN0KQ0KKG1lYW4ocHJlZC5sZGEkY2xhc3MgIT0gbXBnMDEudGVzdCkqMTAwKQ0KYGBgDQogDQpsYSB0YXNhIGRlIGVycm9yIGRlIHBydWViYSBlcyBkZWwgMTIsNiUNCg0KIyMjIyBlKSBSZWFsaWNlIFFEQSBlbiBsb3MgZGF0b3MgZGUgZW50cmVuYW1pZW50byBwYXJhIHByZWRlY2lyICJtcGcwMSIgdXRpbGl6YW5kbyBsYXMgdmFyaWFibGVzIHF1ZSBwYXJlY8OtYW4gZXN0YXIgbcOhcyBhc29jaWFkYXMgY29uICJtcGcwMSIgZW4gKGIpLiDCv0N1w6FsIGVzIGVsIGVycm9yIGRlIHBydWViYSBkZWwgbW9kZWxvIG9idGVuaWRvPw0KDQpgYGB7cn0NCmZpdC5xZGEgPC0gcWRhKG1wZzAxIH4gY3lsaW5kZXJzICsgd2VpZ2h0ICsgZGlzcGxhY2VtZW50ICsgaG9yc2Vwb3dlciwgZGF0YSA9IEF1dG8sDQogICAgICAgICAgICAgICBzdWJzZXQgPSB0cmFpbikNCmZpdC5xZGENCnByZWQucWRhIDwtIHByZWRpY3QoZml0LnFkYSwgQXV0by50ZXN0KQ0KdGFibGUocHJlZC5xZGEkY2xhc3MsIG1wZzAxLnRlc3QpDQoobWVhbihwcmVkLnFkYSRjbGFzcyAhPSBtcGcwMS50ZXN0KSoxMDApDQpgYGANCg0KbGEgdGFzYSBkZSBlcnJvciBkZSBwcnVlYmEgZXMgZGVsIDEzLDE4JQ0KDQojIyMjIGYpIFJlYWxpY2UgdW5hIHJlZ3Jlc2nDs24gbG9nw61zdGljYSBlbiBsb3MgZGF0b3MgZGUgZW50cmVuYW1pZW50byBwYXJhIHByZWRlY2lyICJtcGcwMSIgdXRpbGl6YW5kbyBsYXMgdmFyaWFibGVzIHF1ZSBwYXJlY8OtYW4gZXN0YXIgbcOhcyBhc29jaWFkYXMgY29uICJtcGcwMSIgZW4gKGIpLiDCv0N1w6FsIGVzIGVsIGVycm9yIGRlIHBydWViYSBkZWwgbW9kZWxvIG9idGVuaWRvPw0KDQpgYGB7cn0NCmZpdC5nbG0gPC0gZ2xtKG1wZzAxIH4gY3lsaW5kZXJzICsgd2VpZ2h0ICsgZGlzcGxhY2VtZW50ICsgaG9yc2Vwb3dlciwgZGF0YSA9IEF1dG8sIGZhbWlseSA9IGJpbm9taWFsLCBzdWJzZXQgPSB0cmFpbikNCnN1bW1hcnkoZml0LmdsbSkNCnByb2JzIDwtIHByZWRpY3QoZml0LmdsbSwgQXV0by50ZXN0LCB0eXBlID0gInJlc3BvbnNlIikNCnByZWQuZ2xtIDwtIHJlcCgwLCBsZW5ndGgocHJvYnMpKQ0KcHJlZC5nbG1bcHJvYnMgPiAwLjVdIDwtIDENCnRhYmxlKHByZWQuZ2xtLCBtcGcwMS50ZXN0KQ0KKG1lYW4ocHJlZC5nbG0gIT0gbXBnMDEudGVzdCkqMTAwKQ0KYGBgDQoNCkxhIHRhc2EgZGUgZXJyb3IgZXMgZGVsIDEyLDElDQoNCiMjIyMgZykgRWplY3V0ZSBLTk4gZW4gbG9zIGRhdG9zIGRlIGVudHJlbmFtaWVudG8sIGNvbiB2YXJpb3MgdmFsb3JlcyBkZSBLLCBwYXJhIHByZWRlY2lyIOKAnG1wZzAx4oCdIHV0aWxpemFuZG8gbGFzIHZhcmlhYmxlcyBxdWUgcGFyZWPDrWFuIGVzdGFyIG3DoXMgYXNvY2lhZGFzIGNvbiDigJxtcGcwMeKAnSBlbiAoYikuIMK/UXXDqSBlcnJvcmVzIGRlIHBydWViYSBvYnRpZW5lcz8gwr9RdcOpIHZhbG9yIGRlIEsgcGFyZWNlIGZ1bmNpb25hciBtZWpvciBlbiBlc3RlIGNvbmp1bnRvIGRlIGRhdG9zPw0KDQpgYGB7cn0NCnRyYWluLlggPC0gY2JpbmQoY3lsaW5kZXJzLCB3ZWlnaHQsIGRpc3BsYWNlbWVudCwgaG9yc2Vwb3dlcilbdHJhaW4sIF0NCnRlc3QuWCA8LSBjYmluZChjeWxpbmRlcnMsIHdlaWdodCwgZGlzcGxhY2VtZW50LCBob3JzZXBvd2VyKVshdHJhaW4sIF0NCnRyYWluLm1wZzAxIDwtIG1wZzAxW3RyYWluXQ0Kc2V0LnNlZWQoMSkNCnByZWQua25uIDwtIGtubih0cmFpbi5YLCB0ZXN0LlgsIHRyYWluLm1wZzAxLCBrID0gMSkNCnRhYmxlKHByZWQua25uLCBtcGcwMS50ZXN0KQ0KbWVhbihwcmVkLmtubiAhPSBtcGcwMS50ZXN0KSoxMDANCmBgYA0KDQoNCmBgYHtyfQ0KcHJlZC5rbm4gPC0ga25uKHRyYWluLlgsIHRlc3QuWCwgdHJhaW4ubXBnMDEsIGsgPSAxMCkNCnRhYmxlKHByZWQua25uLCBtcGcwMS50ZXN0KQ0KbWVhbihwcmVkLmtubiAhPSBtcGcwMS50ZXN0KSoxMDANCmBgYA0KDQpgYGB7cn0NCnByZWQua25uIDwtIGtubih0cmFpbi5YLCB0ZXN0LlgsIHRyYWluLm1wZzAxLCBrID0gMTAwKQ0KdGFibGUocHJlZC5rbm4sIG1wZzAxLnRlc3QpDQptZWFuKHByZWQua25uICE9IG1wZzAxLnRlc3QpKjEwMA0KYGBgDQoNCmN1YW5kbyBrPTEwMCBzZSBvYnRpZW5lIGxhIG1lbm9yIHRhc2EgZGUgZXJyb3INCg0KDQojIyMjIDEyKSBFc3RlIHByb2JsZW1hIGltcGxpY2EgZnVuY2lvbmVzIGRlIGVzY3JpdHVyYS4NCg0KIyMjIyBhKSBFc2NyaWJhIHVuYSBmdW5jacOzbiwgUG93ZXIgKCksIHF1ZSBpbXByaW1hIGVsIHJlc3VsdGFkbyBkZSBhdW1lbnRhciAyIGEgbGEgM8KqIHBvdGVuY2lhLiBFbiBvdHJhcyBwYWxhYnJhcywgc3UgZnVuY2nDs24gZGViZSBjYWxjdWxhciAyXjMgZSBpbXByaW1pciBsb3MgcmVzdWx0YWRvcy4NCg0KYGBge3J9DQpQb3dlciA8LSBmdW5jdGlvbigpIHsNCiAgICAyXjMNCn0NCnByaW50KFBvd2VyKCkpDQoNCmBgYA0KDQojIyMjIGIpIENyZWUgdW5hIG51ZXZhIGZ1bmNpw7NuLCBQb3dlcjIoKSwgcXVlIGxlIHBlcm1pdGEgcGFzYXIgY3VhbHF1aWVyYSBkZSBsb3MgZG9zIG7Dum1lcm9zLCAgICJ4IiB5ICJhIiwgZSBpbXByaW1hIGVsIHZhbG9yIGRlICJ4IF4gYSIuDQoNCmBgYHtyfQ0KUG93ZXIyIDwtIGZ1bmN0aW9uKHgsIGEpIHsNCiAgICB4XmENCn0NCg0KUG93ZXIyKDMsIDgpDQpgYGANCg0KIyMjIyBjKSBVc2FuZG8gbGEgZnVuY2nDs24gUG93ZXIyICgpIHF1ZSBhY2FiYSBkZSBlc2NyaWJpciwgY2FsY3VsZSAxMF4zLCA4XjE3IHkgMTMxXjMuDQoNCmBgYHtyfQ0KUG93ZXIyKDEwLCAzKQ0KUG93ZXIyKDgsIDE3KQ0KUG93ZXIyKDEzMSwgMykNCmBgYA0KDQojIyMjIGQpIEFob3JhIGNyZWUgdW5hIG51ZXZhIGZ1bmNpw7NuLCBQb3dlcjMgKCksIHF1ZSByZWFsbWVudGUgZGV2dWVsdmEgZWwgcmVzdWx0YWRvICJ4IF4gYSIgY29tbyB1biBvYmpldG8gUiwgZW4gbHVnYXIgZGUgc2ltcGxlbWVudGUgaW1wcmltaXJsbyBlbiBsYSBwYW50YWxsYS4gRXMgZGVjaXIsIHNpIGFsbWFjZW5hIGVsIHZhbG9yICJ4IF4gYSIgZW4gdW4gb2JqZXRvIGxsYW1hZG8gInJlc3VsdGFkbyIgZGVudHJvIGRlIHN1IGZ1bmNpw7NuLCBlbnRvbmNlcyBzaW1wbGVtZW50ZSBwdWVkZSBkZXZvbHZlciAoKSBlc3RlIHJlc3VsdGFkbw0KDQpgYGB7cn0NClBvd2VyMyA8LSBmdW5jdGlvbih4ICwgYSkgew0KICAgIHJlc3VsdCA8LSB4XmENCiAgICByZXR1cm4ocmVzdWx0KQ0KfQ0KYGBgDQoNCiMjIyMgZSkgQWhvcmEsIHV0aWxpemFuZG8gbGEgZnVuY2nDs24gUG93ZXIzKCksIGNyZWUgdW5hIGdyw6FmaWNhIGRlIGYgKHgpID0geF4zLiBFbCBlamUgeCBkZWJlIG1vc3RyYXIgdW4gcmFuZ28gZGUgZW50ZXJvcyBkZSAxIGEgMTAsIHkgZWwgZWplIHkgZGViZSBtb3N0cmFyIHheMi4gRXRpcXVldGUgbG9zIGVqZXMgYXByb3BpYWRhbWVudGUgeSB1c2UgdW4gdMOtdHVsbyBhcHJvcGlhZG8gcGFyYSBsYSBmaWd1cmEuIENvbnNpZGVyZSBtb3N0cmFyIGVsIGVqZSB4LCBlbCBlamUgeSwgbyBhbWJvcyBlbiBsYSBlc2NhbGEgbG9nYXLDrXRtaWNhLg0KDQpgYGB7cn0NCnggPC0gMToxMA0KcGxvdCh4LCBQb3dlcjMoeCwgMiksIGxvZyA9ICJ4eSIsIHhsYWIgPSAiTG9nIG9mIHgiLCB5bGFiID0gIkxvZyBvZiB4XjIiLCBtYWluID0gIkxvZyBvZiB4XjIgdnMgTG9nIG9mIHgiKQ0KYGBgDQoNCiMjIyMgZikgQ3JlZSB1bmEgZnVuY2nDs24sIFBsb3RQb3dlciAoKSwgcXVlIGxlIHBlcm1pdGEgY3JlYXIgdW4gZ3LDoWZpY28gZGUgIngiIGNvbnRyYSAieCBeIGEiIHBhcmEgdW5hICJhIiBmaWphIHkgcGFyYSB1biByYW5nbyBkZSB2YWxvcmVzIGRlICJ4Ii4NCg0KYGBge3J9DQpQbG90UG93ZXIgPC0gZnVuY3Rpb24oeCwgYSkgew0KICAgIHBsb3QoeCwgUG93ZXIzKHgsIGEpKQ0KfQ0KUGxvdFBvd2VyKDE6MTAsIDMpDQoNCmBgYA0KDQoNCiMjIyMgMTMpIFVzYW5kbyBlbCBjb25qdW50byBkZSBkYXRvcyAiQm9zdG9uIiwgYWp1c3RlIGxvcyBtb2RlbG9zIGRlIGNsYXNpZmljYWNpw7NuIHBhcmEgcHJlZGVjaXIgc2kgdW4gc3VidXJiaW8gZGFkbyB0aWVuZSB1bmEgdGFzYSBkZSBjcmltaW5hbGlkYWQgcG9yIGVuY2ltYSBvIHBvciBkZWJham8gZGUgbGEgbWVkaWFuYS4gRXhwbG9yZSBsb3MgbW9kZWxvcyBkZSByZWdyZXNpw7NuIGxvZ8Otc3RpY2EsIExEQSB5IEtOTiB1dGlsaXphbmRvIHZhcmlvcyBzdWJjb25qdW50b3MgZGUgbG9zIHByZWRpY3RvcmVzLiBEZXNjcmliZSB0dXMgaGFsbGF6Z29zLg0KDQpgYGB7cn0NCmxpYnJhcnkoTUFTUykNCmF0dGFjaChCb3N0b24pDQpjcmltMDEgPC0gcmVwKDAsIGxlbmd0aChjcmltKSkNCmNyaW0wMVtjcmltID4gbWVkaWFuKGNyaW0pXSA8LSAxDQpCb3N0b24gPC0gZGF0YS5mcmFtZShCb3N0b24sIGNyaW0wMSkNCnRyYWluIDwtIDE6KGxlbmd0aChjcmltKSAvIDIpDQp0ZXN0IDwtIChsZW5ndGgoY3JpbSkgLyAyICsgMSk6bGVuZ3RoKGNyaW0pDQpCb3N0b24udHJhaW4gPC0gQm9zdG9uW3RyYWluLCBdDQpCb3N0b24udGVzdCA8LSBCb3N0b25bdGVzdCwgXQ0KY3JpbTAxLnRlc3QgPC0gY3JpbTAxW3Rlc3RdDQpmaXQuZ2xtIDwtIGdsbShjcmltMDEgfiAuIC0gY3JpbTAxIC0gY3JpbSwgZGF0YSA9IEJvc3RvbiwgZmFtaWx5ID0gYmlub21pYWwsIHN1YnNldCA9IHRyYWluKQ0KYGBgDQoNCmBgYHtyfQ0KcHJvYnMgPC0gcHJlZGljdChmaXQuZ2xtLCBCb3N0b24udGVzdCwgdHlwZSA9ICJyZXNwb25zZSIpDQpwcmVkLmdsbSA8LSByZXAoMCwgbGVuZ3RoKHByb2JzKSkNCnByZWQuZ2xtW3Byb2JzID4gMC41XSA8LSAxDQp0YWJsZShwcmVkLmdsbSwgY3JpbTAxLnRlc3QpDQptZWFuKHByZWQuZ2xtICE9IGNyaW0wMS50ZXN0KSoxMDANCmBgYA0KDQpgYGB7cn0NCmZpdC5nbG0gPC0gZ2xtKGNyaW0wMSB+IC4gLSBjcmltMDEgLSBjcmltIC0gY2hhcyAtIG5veCwgZGF0YSA9IEJvc3RvbiwgZmFtaWx5ID0gYmlub21pYWwsIHN1YnNldCA9IHRyYWluKQ0KcHJvYnMgPC0gcHJlZGljdChmaXQuZ2xtLCBCb3N0b24udGVzdCwgdHlwZSA9ICJyZXNwb25zZSIpDQpwcmVkLmdsbSA8LSByZXAoMCwgbGVuZ3RoKHByb2JzKSkNCnByZWQuZ2xtW3Byb2JzID4gMC41XSA8LSAxDQp0YWJsZShwcmVkLmdsbSwgY3JpbTAxLnRlc3QpDQptZWFuKHByZWQuZ2xtICE9IGNyaW0wMS50ZXN0KSoxMDANCg0KYGBgDQoNCmBgYHtyfQ0KZml0LmxkYSA8LSBsZGEoY3JpbTAxIH4gLiAtIGNyaW0wMSAtIGNyaW0sIGRhdGEgPSBCb3N0b24sIHN1YnNldCA9IHRyYWluKQ0KcHJlZC5sZGEgPC0gcHJlZGljdChmaXQubGRhLCBCb3N0b24udGVzdCkNCnRhYmxlKHByZWQubGRhJGNsYXNzLCBjcmltMDEudGVzdCkNCm1lYW4ocHJlZC5sZGEkY2xhc3MgIT0gY3JpbTAxLnRlc3QpKjEwMA0KDQpgYGANCg0KYGBge3J9DQpmaXQubGRhIDwtIGxkYShjcmltMDEgfiAuIC0gY3JpbTAxIC0gY3JpbSAtIGNoYXMgLSBub3gsIGRhdGEgPSBCb3N0b24sIHN1YnNldCA9IHRyYWluKQ0KcHJlZC5sZGEgPC0gcHJlZGljdChmaXQubGRhLCBCb3N0b24udGVzdCkNCnRhYmxlKHByZWQubGRhJGNsYXNzLCBjcmltMDEudGVzdCkNCm1lYW4ocHJlZC5sZGEkY2xhc3MgIT0gY3JpbTAxLnRlc3QpKjEwMA0KDQpgYGANCg0KYGBge3J9DQp0cmFpbi5YIDwtIGNiaW5kKHpuLCBpbmR1cywgY2hhcywgbm94LCBybSwgYWdlLCBkaXMsIHJhZCwgdGF4LCBwdHJhdGlvLCBibGFjaywgbHN0YXQsIG1lZHYpW3RyYWluLCBdDQp0ZXN0LlggPC0gY2JpbmQoem4sIGluZHVzLCBjaGFzLCBub3gsIHJtLCBhZ2UsIGRpcywgcmFkLCB0YXgsIHB0cmF0aW8sIGJsYWNrLCBsc3RhdCwgbWVkdilbdGVzdCwgXQ0KdHJhaW4uY3JpbTAxIDwtIGNyaW0wMVt0cmFpbl0NCnNldC5zZWVkKDEpDQpwcmVkLmtubiA8LSBrbm4odHJhaW4uWCwgdGVzdC5YLCB0cmFpbi5jcmltMDEsIGsgPSAxKQ0KdGFibGUocHJlZC5rbm4sIGNyaW0wMS50ZXN0KQ0KbWVhbihwcmVkLmtubiAhPSBjcmltMDEudGVzdCkqMTAwDQoNCmBgYA0KDQpgYGB7cn0NCnByZWQua25uIDwtIGtubih0cmFpbi5YLCB0ZXN0LlgsIHRyYWluLmNyaW0wMSwgayA9IDEwKQ0KdGFibGUocHJlZC5rbm4sIGNyaW0wMS50ZXN0KQ0KbWVhbihwcmVkLmtubiAhPSBjcmltMDEudGVzdCkqMTAwDQoNCmBgYA0KDQpgYGB7cn0NCnByZWQua25uIDwtIGtubih0cmFpbi5YLCB0ZXN0LlgsIHRyYWluLmNyaW0wMSwgayA9IDEwMCkNCnRhYmxlKHByZWQua25uLCBjcmltMDEudGVzdCkNCm1lYW4ocHJlZC5rbm4gIT0gY3JpbTAxLnRlc3QpKjEwMA0KYGBgDQoNCg0KDQoNCg0KDQo=