III.rmd

---
title: "Métodos y resultados"
author: "Julian Cruz"
csl: apa.csl
output:
  pdf_document:
    keep_tex: yes
    number_sections: yes
    toc: yes
linkcolor: blue
bibliography: Biblio.bib
---

 - Fundamentar el ACP: Por qué y cómo se hace un ACP a las funciones empíricas?

El presente capítulo presenta un estudio en condiciones controladas de las implicaciones prácticas de los componentes teóricos presentados en el capítulo anterior. Siguiendo este orden de ideas es necesario examinar detenidamente las propiedades de las nubes de puntos y de sus posibles semidistancias. Diversas preguntas son respondidas mediante una serie de simlaciones, haciendo uso de la distancia de Minkowski de segundo orden entre funciones, en adelante notada $L_2$; obteniendo como resultado la implementación de los algoritmos J-nubes y K-nubes, que corresponden a la aplicación de conglomerados jerárquicos y de K-means sobre nubes de puntos respectivamente.

El trabajo en condiciones controladas se divide en dos partes, la primera caracteriza la distancia $L_2$ entre funciones, aborda el problema del tamaño y consigue controlarlo generando una semidistancia entre nubes de puntos; la segunda hace uso de esta herramienta para observar el comportamiento de las nubes de puntos al aplicar métodos de agrupación basados en distancias.

Para comodidad de lector la metodología y los resultados son presentadas de manera simultánea.

```{r setup, include=FALSE}
library("reshape2")
library("agricolae")
library("magrittr")
library("Rcpp")
library("FactoMineR")
library("dplyr")
library("ggplot2")
library("fpc")
library("latex2exp")
library("knitr")
library("stargazer")

knitr::opts_chunk$set(echo = FALSE, cache = TRUE)

sourceCpp("distancia.cpp")

```

```{r}
to_nubepuntos <- function(a){
  data.frame(value = a, weight = 1/length(a)) -> nubepuntos
 class(nubepuntos) <- c("data.frame", "nubepuntos")
  nubepuntos
}

```

```{r}
collapsa_nubes <- function(nube_list){
  
  sum(unlist(lapply(nube_list, nrow))) -> n_total
  lapply(nube_list, function(nube){
    nube$weight <- nube$weight*nrow(nube)/n_total
    nube}) -> .
  do.call(rbind, .) -> nubepuntos
    class(nubepuntos) <- c("data.frame", "nubepuntos")
  nubepuntos

}

list_outer <- function(a,b, fun) {
  outer(a, b, function(x,y) vapply(seq_along(x), function(i) fun(x[[i]], y[[i]]), numeric(1)))
}

list_auto <- function(a, fun){list_outer(a, a, fun)}

auto <- function(a, fun){fun(a, a)}

```

# Distancia

El estudio simulado de la distancia entre nubes de puntos parte de la inquietud presentada en el capítulo anterior, en donde el factor de convergencia afecta la construcción presentada de las semidistancias entre nubes de puntos. La semidistancia de Minkowski controlada por tamaño presentada en el capítulo anterior requiere una función $h$ que permita controlar la convergencia a partir del tamaño de las nubes de puntos. Al respecto es necesario examinar esta convergencia y establecer patrones para su comportamiento.

## El problema de tamaño

La convergencia indicada por el teorema de Glivenko Cantelli sucede a medida de que la muestra aumenta. Esto genera efectos sobre la distancia entre nubes de puntos, explicado de manera sucinta con un ejemplo: Dos nubes de puntos de tamaño 10 que tienen una distancia de 1.2 pueden provenir de la misma distribución, sin embargo es muy poco probable que dos nubes de puntos de tamaño 1000 con una distancia igual provengan de la misma distribución. Es posible generar dos nubes pequeñas de la misma distribución y que tengan entre ellas una distancia mayor que dos nubes grandes de distribuciones diferentes como se muestra en la figura HTG.

```{r}
dist_nubes <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
 distan_np(tmp[order(tmp$value),])
}

# rnorm(1000) -> a
# rnorm(840) -> b
# 
# a <- to_nubepuntos(a)
# b <- to_nubepuntos(b)
# 
# dist_nubes(b, a)

```

```{r, fig.height=3, fig.width=3.3}
set.seed(2345)
to_nubepuntos(rnorm(10)) -> n1_peque
to_nubepuntos(rnorm(10)) -> n2_peque

dist_peque <- dist_nubes(n1_peque, n2_peque)

data.frame(n1_peque, n2_peque) %>% ggplot + geom_step(aes(x = value), stat = "ecdf") + geom_step(aes(x = value.1), stat = "ecdf") + theme_minimal() + labs(x = "Valores", y = "FDE")

to_nubepuntos(rnorm(1000, .5)) -> n1_grande
to_nubepuntos(rnorm(1000)) -> n2_grande

dist_grande <- dist_nubes(n1_grande, n2_grande)

data.frame(n1_grande, n2_grande) %>% ggplot + geom_step(aes(x = value), stat = "ecdf") + geom_step(aes(x = value.1), stat = "ecdf") + theme_minimal() + labs(x = "Valores", y = "FDE")

# start.time <- Sys.time()
# dist_nubes(n1_grande, n2_grande)
# end.time <- Sys.time()
# time.taken <- end.time - start.time
# time.taken
# 
# start.time <- Sys.time()
# D_def(n1_grande, n2_grande)
# end.time <- Sys.time()
# time.taken <- end.time - start.time
# time.taken
# 
# start.time <- Sys.time()
# dist_nubes(n1_peque, n2_peque)
# end.time <- Sys.time()
# time.taken <- end.time - start.time
# time.taken
# 
# start.time <- Sys.time()
# D_def(n1_peque, n2_peque)
# end.time <- Sys.time()
# time.taken <- end.time - start.time
# time.taken

```

Figura HTG: Izquerda, funciones empíricas de muestras pequeñas de la misma distribución normal estándar, su distancia es `r round(dist_peque, 2)`; derecha, funciones empíricas de muestras grandes de la distribuciones normales distintas, su distancia es `r round(dist_grande, 2)`.

Si la semidistancia definida refleja características de tamaño en vez de características distribucionales, se trata de un problema de diseño. Por ende es necesario hacer uso de un mecanismo de control que retire la influencia del tamaño en la distancia. Para logar una aproximación adecuada al problema se realiza un examen a cerca de la convergencia de las nubes de puntos.

## Convergencia

El estudio de Monte Carlo respecto a la convergencia de la nube de puntos mide la incidencia del tamaño en la distancia $L^2$ entre nubes de puntos. Se realiza de la siguiente manera:

 - Para la distribución normal estándar.
 - Generar diez mil nubes de puntos: Cien de tamaño 1, cien de tamaño 2, ... cien de tamaño cien.
 - Generar una nube de puntos de tamaño diezmil. Se supone que por su tamaño, esta nube de puntos representa la distribución teórica, o, en otras palabras, una nube de puntos de tamaño infinito.
 - Encontrar la distancia entre cada una de las diezmil nubes de puntos y la nube gigante-teórica.
 - Agrupar estas distancias de acuerdo al tamaño de la nube de puntos que representan y obtener promedios.
 - Graficar las distancias meias por tamaño y ajustar una curva opcional.
 - Hacer lo mismo para la distribución uniforme.

```{r, fig.height = 3, fig.width = 3.3}
lapply(1:100, function(x){
  1:100 %>% lapply(function(y){rnorm(n = x, mean = 0, sd = 1)}) %>% lapply(to_nubepuntos)
  }) -> nubes1

nube_grande1 <- to_nubepuntos(rnorm(10000))

nubes1 %>% lapply(list_outer, b = list(nube_grande1), dist_nubes) %>% lapply(mean) %>% unlist -> distancias1

lm(log(distancias1) ~ log(1:100)) %$% coefficients -> coef1

distancias1 %>% 
  data.frame(X = 1:100, Y = exp(coef1[1])*(1:100)^coef1[2]) %>% 
  ggplot + geom_point(aes(x = X, y = .)) +
  geom_line(aes(x = X, y = Y), colour = "#0000aa") + theme_minimal() + 
  labs(x = "Tamaño", y = "Distancia media")

#, caption = TeX("Distancias medias por tamaño para\nuna distribución normal estándar: $h(x) \\propto \\frac{1}{n_s}$")

lapply(1:100, function(x){
  1:100 %>% lapply(function(y){runif(n = x)}) %>% lapply(to_nubepuntos)
  }) -> nubes2

nube_grande2 <- to_nubepuntos(runif(10000))

nubes2 %>% lapply(list_outer, b = list(nube_grande2), dist_nubes) %>% lapply(mean) %>% unlist -> distancias2

lm(log(distancias2) ~ log(1:100)) %$% coefficients -> coef2

distancias2 %>% 
  data.frame(X = 1:100, Y = exp(coef2[1])*(1:100)^coef2[2]) %>% 
  ggplot + geom_point(aes(x = X, y = .)) + 
  geom_line(aes(x = X, y = Y), colour = "#0000aa") + theme_minimal() + 
  labs(x = "Tamaño", y = "Distancia media")
# , caption = TeX("Distancias medias por tamaño para\nuna distribución uniforme: $h(x) \\propto \\frac{1}{n_s}$")
```

Figura QRG: Izquierda, Distancias medias por tamaño para una distribución normal; Distancias medias por tamaño para una distribución uniforme. Una primera aproximación intuitiva se obtiene al ajustar una curva inversamente proporcional al tamaño de la nube de puntos.

A partir de los resultados es posible observar como la distancia media por tamaño sigue un ajuste inversamente proporcional al tamaño de la nube de puntos (Figura QRG). Este procedimiento da pautas respecto a la mejor función para corregir la distancia $L_2$.

## Corrección por tamaño

La tarea de corregir la distancia por tamaño consiste en encontrar una función $h$ que al multiplicarla por la distancia $L^2$ retire total o parcialmente la incidencia del tamaño de las nubes de puntos en la semidistancia propuesta. Al respecto las funciones candidatas son: $h(n_{s_1}, n_{s_2}) = 1$, $h(n_{s_1}, n_{s_2}) = n_{s_1} + n_{s_2}$, $h(n_{s_1}, n_{s_2}) = \max(n_{s_1}, n_{s_2})$, $h(n_{s_1}, n_{s_2}) = (n_{s_1} + n_{s_2})^2$, $h(n_{s_1}, n_{s_2}) = \min(n_{s_1}, n_{s_2})$ y $h(n_{s_1}, n_{s_2}) = \sqrt{\min(n_{s_1}, n_{s_2})}$.

La simulación se realiza de la manera siguiente:

```{r}
sample_size <- 100

```

 - Generar `r sample_size` nubes de puntos de distribución normal (uniforme).
 - Generar la matriz de distancias para cada una de las correcciones $h$ mencionadas.
 - Se muestra la matriz de distancias.
 - Se busca aquella que sea más homogénea.

```{r}
D_EMD <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
  EMD(tmp[order(tmp$value),])
}

D0 <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
  distan_np(tmp[order(tmp$value),])
}

D1 <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
  nrow(tmp) * distan_np(tmp[order(tmp$value),])
}

D2 <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
  max(nrow(nube_0), nrow(nube_1)) * distan_np(tmp[order(tmp$value),])
}

D3 <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
  nrow(tmp) * nrow(tmp) * distan_np(tmp[order(tmp$value),])
}

D4 <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
  min(nrow(nube_0), nrow(nube_1)) * distan_np(tmp[order(tmp$value),])
}

D5 <- function(nube_0, nube_1){
  rbind(
    cbind(nube_0, binary = 0),
    cbind(nube_1, binary = 1)
  ) -> tmp
  
  sqrt(min(nrow(nube_0), nrow(nube_1))) * distan_np(tmp[order(tmp$value),])
}

```

```{r}
lapply(1:sample_size, rnorm) %>% lapply(to_nubepuntos) -> nubes_norm
lapply(1:sample_size, runif) %>% lapply(to_nubepuntos) -> nubes_unif

lapply(
  list(D0 = D0, D1 = D1, D2 = D2, D3 = D3, D4 = D4, D5 = D5),
  function(fun){
    lapply(
      list(normal = nubes_norm, uniforme = nubes_unif),
      list_auto,
      fun = fun
    )
  }
  ) %>% do.call(c, .) %>% lapply(as.vector) %>% data.frame(expand.grid(1:sample_size,1:sample_size)) -> dis_frame

```

```{r, fig.width = 3.5, fig.height = 2.5}
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D0.normal) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D0.uniforme) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")

```

Figura HRR: Distancia corregida mediante $h(n_{s_1}, n_{s_2}) = 1$. Esta es la visualización de las matrices de distancias sin controlar el tamaño. Como se mencionó, los tamaños de muestra pequeños producen distancias mayores a pesar de provenir de la misma distribución. A la izquierda, se usa una distribución normal estándar; a la derecha, una distribución uniforme.


```{r, fig.width = 3.5, fig.height = 2.5}
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D1.normal) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D1.uniforme) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")

```

Figura HRJ: Distancia corregida mediante $h(n_{s_1}, n_{s_2}) = n_{s_1} + n_{s_2}$. Visualización de la matriz de distancias corrigiendo por la suma de los tamaños. Es notorio como las distancias crecen cuando los tamaños de las nubes de puntos son muy distintos a pesar de provenir de la misma distribución. A la izquierda, se usa una distribución normal estándar; a la derecha, una distribución uniforme.

```{r, fig.width = 3.5, fig.height = 2.5}
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D2.normal) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D2.uniforme) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")

```

Figura HRV: Distancia corregida mediante $h(n_{s_1}, n_{s_2}) = \max(n_{s_1}, n_{s_2})$. Visualización de la matriz de distancias corrigiendo por el máximo entre los tamaños. Es notorio como las distancias crecen cuando los tamaños de las nubes de puntos son muy distintos a pesar de provenir de la misma distribución. A la izquierda, se usa una distribución normal estándar; a la derecha, una distribución uniforme.

```{r, fig.width = 3.5, fig.height = 2.5}
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D3.normal) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D3.uniforme) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")

```

Figura HRQ: Distancia corregida mediante $h(n_{s_1}, n_{s_2}) = (n_{s_1} + n_{s_2})^2$. Visualización de la matriz de distancias corregida la suma de los tamaños al cuadrado. A la izquierda, se usa una distribución normal estándar; a la derecha, una distribución uniforme.

```{r, fig.width = 3.5, fig.height = 2.5}
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D4.normal) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D4.uniforme) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia")+ scale_fill_gradient(high = "#333355", low = "#ffffff")

```

Figura HRL: Distancia corregida mediante $h(n_{s_1}, n_{s_2}) = \min(n_{s_1}, n_{s_2})$. Visualización de la matriz de distancias corregida mediante el mínimo de los tamaños. Los valores altos, que se encontraban en un principio asociados a tamaños bajos ahora se relacionan con tamaños altos. A la izquierda, se usa una distribución normal estándar; a la derecha, una distribución uniforme.

```{r, fig.width = 3.5, fig.height = 2.5}
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D5.normal) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia") + scale_fill_gradient(high = "#333355", low = "#ffffff")
dis_frame %>% ggplot + aes(x = Var1, y = Var2, fill = D5.uniforme) + geom_raster() + theme_minimal(base_size = 6) + labs(x ="", y = "", fill = "Distancia") + scale_fill_gradient(high = "#333355", low = "#ffffff")

```

Figura HRJ: Distancia corregida mediante $h(n_{s_1}, n_{s_2}) = \sqrt{\min(n_{s_1}, n_{s_2})}$. Visualización de la matriz de distancias corrigiendo mediante la raíz cuadrada del mínimo de los tamaños. La distribución aleatoria de valores altos y bajos a través de todos los tamaños es notoria. A la izquierda, se usa una distribución normal estándar; a la derecha, una distribución uniforme.

El análisis de las visualizaciones anteriores conduce a pensar que la corrección usando la raíz cuadrada del mínimo de los tamaños es la más eficiente para estabilizar la distancia en torno a la convergencia, resolviendo el problema planteado. Esto es consistente particularmente con la función ajustada en la sección anterior, donde se examinaba la distancia media según el tamaño. Para corroborar este hallazgo se plantea el modelo de regresión:

$$ln(d(s_1, s_2)) = \beta_0 + \beta_1ln(\min(n_{s_1}, n_{s_2})) + \beta_2ln(\max(n_{s_1}, n_{s_2}))$$
Que se presenta para datos normales:

```{r}
dis_frame %>% filter(Var1 != Var2) %>% mutate(min_12 = apply(cbind(Var1, Var2), 1, min), max_12 = apply(cbind(Var1, Var2), 1, max)) %>% lm(log(D0.normal) ~ log(min_12) + log(max_12), data = .) -> modelo_1

stargazer(modelo_1, summary = TRUE)

modelo_1 %>% summary

```

Y para datos uniformes:

```{r}
dis_frame %>% filter(Var1 != Var2) %>% mutate(min_12 = apply(cbind(Var1, Var2), 1, min), max_12 = apply(cbind(Var1, Var2), 1, max)) %>% lm(log(D0.uniforme) ~ log(min_12) + log(max_12), data = .) -> modelo_2

stargazer(modelo_2, summary = TRUE)

modelo_2 %>% summary

```

Ambos modelos ratifican la cercanía de $\beta_1$ a $\frac{-1}{2}$ además de la posibilidad de pasar por alto una corrección con respecto al máximo. En este sentido la distancia a usar queda definida así:^[Bueno, acá es donde toca explicar la cosa]

## Criterio de Cramer-von Misses generalizado controlado por tamaño entre nubes de puntos

Sean $(\Omega, \mathcal{B})$. Dadas dos nubes de puntos $s_1$ y $s_2$, con sus medidas de probabilidad empírica $P_{s_1}$ y $P_{s_2}$ y sus funciones de densidad general empírica $f_{s_1}$ y $f_{s_2}$ y tamaños $n_{s_1}$ y $n_{s_2}$; se define la función 

$$\mathtt{d}(s_1, s_2) =  \min\left(\sqrt{\min(n_{s_1}, n_{s_2})\int(f_{s_1} - f_{s_2})^2 d\lambda}, \sqrt{\int(f_{s_1} - f_{s_2})^2 d\lambda} + 1\right)$$

Como ya se mencionó, $\mathtt{d}$ no es una distancia: nubes de puntos equivalentes tienen medidas de pribabilidad empírica iguales y funciones de densidad general empírica iguales, por lo cual dos nubes de puntos distintas pueden tener distancia cero. Tampoco es una semidistancia, pues $\sqrt{\min(n_{s_1}, n_{s_2})}$ no es constante, en conescuencia es posible que no cumpla con la desigualdad triangular en unos casos muy específicos. El caso contínuo puede examinarse a la luz de las funciones de densidad acumulada y acumulada empírica, no obstante los resultados teóricos son bastante más generales.

# Métodos de agrupación

La estructura de similitud entre objetos permite la aplicación de algoritmos de agrupación basados en distancias, semidistancias y criterios de similitud. Haciendo uso del criterio definido se estructuran algoritmos de agrupación entre nubes de puntos: J-nubes, que hace referencia a algoritmos de agrupación jerárquica y K-nubes, basado en K-means. Ambos algoritmos deberían agrupar las nubes de puntos de acuerdo a similitudes en su medida de probabilidad empírica, reflejando estructuras distribucionales más que características paramétricas.

En esta sección se presenta un estudio de Monte Carlo comparativo que muestra resultados de los métodos propuestos J-nubes y K-nubes en paralelo con abordajes más tradicionales como las prueba de Duncan y LSD de Fisher. Para tal fin se simulan 200 datos para 20 tratamientos con distribución $N(\mu_i, \sigma_i)$ según la tabla~\ref{T3}.

\begin{table}[htb!]
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline 
Nube & $\mu$ & $\sigma$ & Nube & $\mu$ & $\sigma$ \\ 
\hline 
A & 1 & $S$ & K & 3 & $S$ \\ 
\hline 
B & 1 & $S$ & L & 3 & $S$ \\ 
\hline 
C & 1 & $S$ & M & 4 & $5S$ \\ 
\hline 
D & 1 & $S$ & N & 4 & $5S$ \\ 
\hline 
E & 2 & $3S$ & O & 4 & $5S$ \\ 
\hline 
F & 2 & $3S$ & P & 4 & $5S$ \\ 
\hline 
G & 2 & $3S$ & Q & 5 & $S$ \\ 
\hline 
H & 2 & $3S$ & R & 5 & $S$ \\ 
\hline 
I & 3 & $S$ & S & 5 & $S$ \\ 
\hline 
J & 3 & $S$ & T & 5 & $S$ \\ 
\hline 
\end{tabular} 
\caption{Simulación de 20 nubes de puntos a partir de una distribución $N(\mu_i,\sigma_i)$, cada uno con 10 datos.}
\label{T3}
\end{table}

```{r}
S <- seq(2, 10, by = 2)

```

Se trata de 20 nubes de puntos con distribución normal con distintas medias y varianzas. La constante $S$, permite controlar la dispersión de los datos, así al variar $S$ se controla la dificultad para distinguir entre las distintas nubes de puntos; de esta manera son realizados `r length(S)` experimentos con $S$ en $\{`r S`\}$. Los métodos presentados deben clasificar las nubes de puntos tal como se simularon o que por lo menos acercarse. Las pruebas de Duncan y LSD de Fisher son usadas para tener un panorama comparativo. En todos los casos se eligió una clasificación en $K=5$ clases.

A continuación se incluyen los resultados para $S = `r S[1]`$ y $S = `r S[5]`$, los resultados restantes se incluyen al final del capítulo.

```{r}
ep_means <- function(lista_nubes, n_grupos = NULL, grupos_iniciales = NULL){
  
  if(is.null(grupos_iniciales)) {
    sample(length(lista_nubes), n_grupos, replace = FALSE) -> cuales
    
    centros <- setNames(lista_nubes[cuales], LETTERS[1:n_grupos])
  } else{
    split(lista_nubes, grupos_iniciales) -> agrupados
    lapply(agrupados, collapsa_nubes) -> centros
    
  }
  
  antes <- 0
  ahora <- 1000
  umbral <- 0.0002
  while(abs(antes - ahora) > umbral){
    antes <- ahora
    list_outer(centros, lista_nubes, D_EMD) -> dist_matriz
    letters[apply(dist_matriz, 2, which.min)] -> grupos
    sum(apply(dist_matriz, 2, min)) -> ahora
    split(lista_nubes, grupos) -> agrupados
    lapply(agrupados, collapsa_nubes) -> centros
    
  }
  
  data.frame(trt = names(lista_nubes), Grupo = grupos, stringsAsFactors = FALSE)
  
}


```

```{r}
k_clouds <- function(lista_nubes, n_grupos = NULL, grupos_iniciales = NULL){
  
  if(is.null(grupos_iniciales)) {
    sample(length(lista_nubes), n_grupos, replace = FALSE) -> cuales
    
    centros <- setNames(lista_nubes[cuales], LETTERS[1:n_grupos])
  } else{
    split(lista_nubes, grupos_iniciales) -> agrupados
    lapply(agrupados, collapsa_nubes) -> centros
    
  }
  
  antes <- 0
  ahora <- 1000
  umbral <- 0.0002
  while(abs(antes - ahora) > umbral){
    antes <- ahora
    list_outer(centros, lista_nubes, D5) -> dist_matriz
    letters[apply(dist_matriz, 2, which.min)] -> grupos
    sum(apply(dist_matriz, 2, min)) -> ahora
    split(lista_nubes, grupos) -> agrupados
    lapply(agrupados, collapsa_nubes) -> centros
    
  }
  
  data.frame(trt = names(lista_nubes), Grupo = grupos, stringsAsFactors = FALSE)
  
}

# ep_means(lista_nubes, 5)
# 
# lapply(1:100, function(i){rnorm(i)}) %>% lapply(to_nubepuntos) -> datos
# 
# lapply(1:10, function(i){rnorm(rpois(1, 100), 50)}) %>% lapply(to_nubepuntos) -> centros
# 
# list_outer(datos, datos, dist_nubes) -> ppp

```

```{r}
evaluated_ecdf <- function(categories, values){

  # a vector with sorted and unique values
  unique(values) -> .
  sort(.) -> ordenados

  # a list of functions, one function per category
  ecdf_categories <- tapply(X = values, INDEX = categories, FUN = ecdf)

  # a data.frame with the evaluations of the functions in the sorted values
  sapply(X = ecdf_categories, FUN = function(x, y){x(y)},
    y = ordenados) -> .
  t(.) ->.
  data.frame(.)->.
  set_colnames(., ordenados)
}

```


```{r}
# setwd("/media/TI106180W0A/Users/julian cruz/Desktop/Investigaciones/CLASIFICA/R")
# source("CLASS.R")

Simul_tabla <- function(sigma){
  
  data.frame(
    Y = c(
      rnorm(40,10,sigma),
      rnorm(40,20,3*sigma),
      rnorm(40, 30,sigma), 
      rnorm(40, 40, 5*sigma), 
      rnorm(40, 50, sigma)
    ),
    X = rep(
      LETTERS[1:20], 
      each=10),
    grupo = rep(letters[1:5], each = 40),
    stringsAsFactors = FALSE
  ) -> datos
  
  datos %>% ggplot + aes(x = Y, colour = grupo, group = X) + geom_step(stat = "ecdf") + theme_minimal() -> gr
  
  # Media
  Media_trt <- datos %$% aggregate(Y, list(X), mean) %>% setNames(c("trt", "Media"))
  
  #Duncan
  aov(Y~X, data = datos) %>% duncan.test("X") %$% groups %>% mutate(trt = rownames(.), grupos = groups) %>% select(-Y, -groups) -> DUNCAN
  
  #LSD
  aov(Y~X, data = datos) %>% LSD.test("X") %$% groups %>% mutate(trt = rownames(.)) %>% select(-Y) -> LSD
  
  #K nubes
  datos %$% split(Y, X) %>% lapply(to_nubepuntos) %>% k_clouds(n_grupos = 5) -> CC

  #J nubes
  datos %$% split(Y, X) %>% lapply(to_nubepuntos) %>% list_auto(D5) %>% as.dist %>% hclust(method = "ward.D") %>% cutree(5) %>% "["(letters, .) -> CC$Ward

  #EP Means
  datos %$% split(Y, X) %>% lapply(to_nubepuntos) %>% ep_means(n_grupos = 5) %$% Grupo -> CC$EP


  list(
    Media_trt,
    data.frame(
      trt = LETTERS[1:20],
      Grupo = rep(letters[1:5], each = 4),
      stringsAsFactors = FALSE
    ),
    CC,
    # LSD$groups[-2], 
    # DUNCAN$groups[-2]
    LSD,
    DUNCAN
  ) -> tmp
  
  Reduce(function(dtf1,dtf2) left_join(dtf1,dtf2, by = "trt"), tmp) %>% 
    setNames(c("Tratamiento", "Media", "Grupo Inicial", "K Nubes", "J Nubes", "EP Means", "LSD", "Duncan")) -> tabla
  
datos %$% evaluated_ecdf(X, Y) -> ecdf_datos

ecdf_datos %>% PCA(graph = FALSE) %$% ind %$% coord %>% as.data.frame %>% cbind(unique(datos[c("X", "grupo")])) %>% ggplot + aes(x = Dim.1, y = Dim.2, colour = grupo) + geom_point()  + theme_minimal() -> pca_gr

  list(gr, tabla, pca_gr, ecdf_datos)
  
}

# sigma <- 2

```


```{r}
set.seed(34567)

lapply(S, Simul_tabla) -> Resultados

```

Para $S = `r S[1]`$ la variabilidad de los datos es pequeña con respecto a la variabilidad total. Los grupos de funciones empíricas pueden distinguirse fácilmente.

```{r, fig.width = 8, fig.height = 2.5}
Resultados[[1]][[1]]

```

El resultado de los distintos métodos hace patente la diferencia entre los grupos. La baja variabilidad de las nubes de puntos facilita la tarea de clasificación.

```{r}
kable(Resultados[[1]][[2]])

```

La visualización mediante componentes principales evidencia los grupos y facilita la interpretación.

```{r, fig.width = 4, fig.height = 2.5}
Resultados[[1]][[3]]

```

Para $S = `r S[5]`$ la variabilidad de los datos es alta con respecto a la variabilidad total, de manera que no resulta facil distinguir a simple vista las funciones de densidad acumulada empíricas coorrespondientes a los distintos grupos de nubes de puntos. Asímismo los métodos empleados tienen inconvenientes en la detección de dichos grupos.

```{r, fig.width = 8, fig.height = 2.5}
Resultados[[5]][[1]]

```

Las dificultades en la clasificacción son notorias, es posible observar un desempeño de J-nubes bastante mejor que el de K-nubes, las pruebas de Duncan y LSD de Fisher presentan grandes dificultades, en éstas, las nubes de puntos con letras en común no presentan diferencias estadísticamente significativas.

```{r}
kable(Resultados[[5]][[2]])

```

La visualización mediante componentes principales difícilmente evidencia los grupos.

```{r, fig.width = 4, fig.height = 2.5}
Resultados[[5]][[3]]

```

Como se ve en as simulaciones realizadas, el procedimiento que mejor clasifica las nubes de puntos es J Nubes. No obstante es posible que realizar una combinación de os algoritmos H Nubes y J Nubes desemboque en un mejor esempeño. En los siguientes capítulos se aplican estos algoritmos a datos reales, estableciendo una utilidad práctica para la agrupación de nubes de puntos.

# Anexo: Resultados

```{r}
para_ecsago <- paste("simulacion_", 1:5, ".csv", sep = "")

i <- 1
```

## $S = `r S[i]`$

```{r, fig.width = 8, fig.height = 2.5}
Resultados[[i]][[1]]

```

```{r}
kable(Resultados[[i]][[2]])

```

```{r, fig.width = 4, fig.height = 2.5}
Resultados[[i]][[3]]

write.table(Resultados[[i]][[4]], para_ecsago[i], sep = ";", row.names = FALSE, dec = ",")

```

```{r}
i <- 2
```

## $S = `r S[i]`$

```{r, fig.width = 8, fig.height = 2.5}
Resultados[[i]][[1]]

```

```{r}
kable(Resultados[[i]][[2]])

```

```{r, fig.width = 4, fig.height = 2.5}
Resultados[[i]][[3]]

write.table(Resultados[[i]][[4]], para_ecsago[i], sep = ";", row.names = FALSE, dec = ",")

```

```{r}
i <- 3
```

## $S = `r S[i]`$

```{r, fig.width = 8, fig.height = 2.5}
Resultados[[i]][[1]]

```

```{r}
kable(Resultados[[i]][[2]])

```

```{r, fig.width = 4, fig.height = 2.5}
Resultados[[i]][[3]]

write.table(Resultados[[i]][[4]], para_ecsago[i], sep = ";", row.names = FALSE, dec = ",")

```

```{r}
i <- 4
```

## $S = `r S[i]`$

```{r, fig.width = 8, fig.height = 2.5}
Resultados[[i]][[1]]

```

```{r}
kable(Resultados[[i]][[2]])

```

```{r, fig.width = 4, fig.height = 2.5}
Resultados[[i]][[3]]

write.table(Resultados[[i]][[4]], para_ecsago[i], sep = ";", row.names = FALSE, dec = ",")

```

```{r}
i <- 5
```

## $S = `r S[i]`$

```{r, fig.width = 8, fig.height = 2.5}
Resultados[[i]][[1]]

```

```{r}
kable(Resultados[[i]][[2]])

```

```{r, fig.width = 4, fig.height = 2.5}
Resultados[[i]][[3]]

write.table(Resultados[[i]][[4]], para_ecsago[i], sep = ";", row.names = FALSE, dec = ",")

```