En este documento se describen algunas de las funcionalidades del paquete dirichletprocess, el cual permite realizar análisis bayesiano no paramétrico haciendo uso de procesos Dirichlet, lo que implica que los usuarios de R no tienen que crear sus propios algoritmos sino que pueden utilizar estas funcionalidades genéricas.

Los autores de este reciente paquete en R son Gordon J. Ross y Dean Markwick. Para información más detallada visitar https://cran.r-project.org/web/packages/dirichletprocess/vignettes/dirichletprocess.pdf.

Instalación y carga del paquete

La instalación y carga del paquete dirichletprocess se realiza de la manera habitual.

# install.packages("dirichletprocess")
library(dirichletprocess)

Estimación de densidades

Ejemplo 1

Se desea realizar la estimación de la densidad de un conjunto de datos simulados:

set.seed(123) # Establecimiento de una semilla
y <- rt(200, 3) + 2 # Generación de los datos muestrales a través de la distribución T de Student
head(y) # Observación de algunos de los datos simulados
## [1]  1.2242842  0.7224712  1.8582340 -0.3597923  3.6535141  4.4295408

Para los datos almacenados en la variable “y”, se realiza la estimación de la densidad de los datos haciendo uso de un modelo de mezcla del proceso Dirichlet de distribuciones Gaussianas:

dp <- DirichletProcessGaussian(y) # Función propia del paquete estudiado que permite crear el modelo de mezcla de Gaussianas univariadas.

dp # Resultados más relevantes de la función anterior
## Dirichlet process object run for 0 iterations.
##                                      
##   Mixing distribution          normal
##   Base measure parameters  0, 1, 1, 1
##   Alpha Prior parameters         2, 4
##   Conjugacy                 conjugate
##   Sample size                     200
##

Posterior al establecimiento de la kernel del modelo de mezcla (Gaussiana con media y varianza desconocidas) y a los parámetros de la medida base del proceso Dirichlet (en este caso los valores por defecto), se hace uso de la función Fit para inferir sobre los parámetros de cada clúster:

dp <- Fit(dp, 1000) # Se realizan 1000 iteraciones para el objeto dp

dp # Resultados actualizados para el objeto dp
## Dirichlet process object run for 1000 iterations.
##                                      
##   Mixing distribution          normal
##   Base measure parameters  0, 1, 1, 1
##   Alpha Prior parameters         2, 4
##   Conjugacy                 conjugate
##   Sample size                     200
##                                      
##   Mean number of clusters        9.63
##   Median alpha                   1.16

El objeto “dp”, al ser ajustado por la función Fit, devuelve ciertas muestras posteriores que pueden ser utilizadas para realizar inferencia:

dp$clusterParametersChain[[10]][[1]] # Muestras MCMC de los parámetros de clúster
## , , 1
## 
##          [,1]
## [1,] 2.086693
## 
## , , 2
## 
##          [,1]
## [1,] 3.233543
## 
## , , 3
## 
##           [,1]
## [1,] -0.589585
## 
## , , 4
## 
##             [,1]
## [1,] -0.05734715
## 
## , , 5
## 
##           [,1]
## [1,] -2.867592
## 
## , , 6
## 
##           [,1]
## [1,] 0.5485425
## 
## , , 7
## 
##            [,1]
## [1,] -0.3851703
## 
## , , 8
## 
##           [,1]
## [1,] 0.6890181
## 
## , , 9
## 
##           [,1]
## [1,] 0.3488289
dp$weightsChain[[10]] # Muestras MCMC de los pesos asociados
## [1] 0.805 0.100 0.030 0.025 0.010 0.010 0.010 0.005 0.005
dp$alphaChain[1:10] # Muestras MCMC del parámetro de concentración
##  [1] 0.1615142 0.2464855 0.5699669 0.3970536 0.8810918 1.1160247 1.0852115
##  [8] 1.7152331 1.7066049 1.4940442

Ejemplo 2

Para este ejemplo se hará uso del conjunto de datos de R denominado Faithful, el cual presenta el tiempo de erupción (en minutos) y los tiempos de espera entre erupciones (en minutos) del géiser Old Faithful en el Parque Nacional de Yellowstone en Wyoming, Estados Unidos. Nos enfocaremos únicamente en la estimación de la densidad de los tiempos de espera.

head(faithful) # Muestra del conjunto de datos
##   eruptions waiting
## 1     3.600      79
## 2     1.800      54
## 3     3.333      74
## 4     2.283      62
## 5     4.533      85
## 6     2.883      55
faithfulTransformed <- scale(faithful$waiting) # Se transforman los tiempos de espera para que tengan media cero y desviación estándar uno

Ya con el conjunto de datos listo, se modelan los tiempos de espera como una mezcla de distribuciones Gaussianas.

its <- 500 # Número de iteraciones
dp <- DirichletProcessGaussian(faithfulTransformed) # Modelo de mezclas del proceso Dirichlet con una kernel Gaussiana univariada
dp <- Fit(dp, its) # Ajuste del modelo con 500 iteraciones
plot(dp) # Gráfico de la densidad de los datos y la estimación de la densidad a partir de la media posteriori del MMPD

A partir del gráfico anterior se puede observar que la naturaleza multi-modal de la densidad original de los datos (en color negro) es capturada de buena manera a través del modelo de mezcla del proceso Dirichlet (MMPD) utilizando una kernel Gaussiana. La línea continua en color rojo representa la estimación de la densidad a través de dicho modelo, mientras que los intervalos de credibilidad se encuentran superpuestos por las líneas punteadas.

Ahora bien, la configuración por defecto del argumento data_method de la función plot es de estimación de la densidad. En ocasiones, la cantidad de datos es bastante limitada, por lo que los autores recomiendan el uso de un histograma.

plot(dp, data_method="hist")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Estimación de la densidad en intervalos acotados

Para este ejemplo se plantea la situación donde se quiera estimar densidades en intervalos delimitados, es decir, suponiendo que se tiene conocimiento de que las \(y_i\) se mueven por ejemplo en el intervalo [0,1].

Afirmar que es posible hacer uso de las mezclas de distribuciones gaussianas sería algo inapropiado pero sí se podría pensar en hacer uso de mezclas de distribuciones Beta, pues se presenta masa positiva en dicho intervalo.

Ahora, se procede a generar la mezcla de dos distribuciones Beta:

set.seed(1000) # Establecimiento de la semilla
y <- c(rbeta(150, 1, 3), rbeta(150, 7, 3)) #Generación de la muestra de datos de las dos distribuciones Beta.
dp <- DirichletProcessBeta(y, 1) #Modelo de mezclas del proceso Dirichlet con una kernel beta bivariada
## Accept Ratio:  0.001
iteraciones=1000 #Número de iteraciones a realizar
dp <- Fit(dp, iteraciones) #Ajuste del modelo con 1000 iteraciones.

Ahora, como se quiere comparar la distribución posteriori resultante con una distribución conocida, se deben realizar las simulaciones manualmente de las muestras posteriores, calcular los intervalos de credibilidad y graficar los resultados.

posteriorFrame <- PosteriorFrame(dp, ppoints(100), ci_size = 0.05) # Calcula la media posteriori y los cuantiles del modelo de mezclas Dp. Simulación.
trueFrame <- data.frame(x=ppoints(100), y=0.5*dbeta(ppoints(100), 1,3)+0.5*dbeta(ppoints(100), 7, 3)) # dataframe donde x representa la secuencia de puntos de probabilidad y Y representa la distribución posteriori de la mezcla de las dos distribuciones beta. Datos verdaderos.

library(ggplot2)
#Gráfico para las densidades
ggplot() +
    geom_ribbon(data=posteriorFrame, #Intervalos de credibilidad para densidad obtenida por simulación.
                  aes(x=x, ymin=X2.5., ymax=X97.5.),
                  alpha=0.2,
                  colour=NA,
                  fill="red") +
    geom_line(data=posteriorFrame, aes(x=x, y=Mean), colour="red") + #Densidad posteriori obtenida por las funciones del paquete. Simulación.
    geom_line(data=trueFrame, aes(x=x, y=y)) #Densidad verdadera.

En la gráfica anterior, la curva negra representa la densidad para la media verdadera posteriori, la roja representa la simulación de la media posteriori y las bandas rosadas representa los intervalos de credibilidad, es posible observar que las dos curvas presentan un comportamiento similar y por tanto, la naturaleza de la densidad original es bien capturada a través del modelo de mezcla del proceso Dirichlet utilizando una kernel Beta Bivariada.

Análisis (multivariado) de clústeres

En concordancia con lo expuesto por los autores del paquete, los procesos Dirichlet proveen una forma efectiva de realizar análisis de clúster, puesto que cada punto muestral \(y_{i}\) es asignado a un parámetro único de clúster \(\theta_{i}\), lo cual es una forma evidente de agrupación de datos.

Para este ejemplo se hará uso de las etiquetas de clúster (clusterLabels dentro del paquete), para indicar a qué grupo pertenece cada uno de los datos de la muestra. Para mostrar esta funcionalidad del paquete se hará uso nuevamente del conjunto de datos Faithful:

faithfulTrans <- scale(faithful) # Se transforman los tiempos de espera para que tengan media cero y desviación estándar uno

head(faithfulTrans) # Muestra de los datos transformados
##     eruptions    waiting
## 1  0.09831763  0.5960248
## 2 -1.47873278 -1.2428901
## 3 -0.13561152  0.2282418
## 4 -1.05555759 -0.6544374
## 5  0.91575542  1.0373644
## 6 -0.52987412 -1.1693335

Como se tienen datos de dos dimensiones, se establecerá un objeto que tenga en consideración una mezcla de distribuciones Gaussianas multivariantes:

set.seed(123) # Establecimiento de la semilla
dp4 <- DirichletProcessMvnormal(faithfulTrans) # Mezcla Dirichlet de distribuciones Gaussianas multivariadas

dp4 <- Fit(dp4, 1000) # 1000 iteraciones a través del algoritmo de Cadenas de Markov de Monte Carlo (MCMC)

Del objeto “dp4” es posible extraer las etiquetas de clúster, las cuales indicarán el grupo al que pertenece cada uno de los datos:

dp4$clusterLabels # Etiquetas de clúster para cada par de datos del conjunto de datos Faithful
##   [1] 2 1 2 1 2 1 2 2 1 2 1 2 2 1 2 1 1 2 1 2 1 1 2 1 2 2 1 2 2 2 2 2 2 2 2 1 1
##  [38] 2 1 2 2 1 2 1 2 2 2 1 2 1 2 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 2 2 1 2 2 1 2 2
##  [75] 1 2 1 2 2 2 2 2 2 1 2 2 2 2 1 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2 2 1 2 1 2 2 2
## [112] 1 2 2 1 2 1 2 1 2 1 2 2 1 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 2 2 1 2 1
## [149] 2 1 2 2 1 2 2 2 2 2 1 2 1 2 1 2 2 2 1 2 1 2 1 1 2 2 2 2 2 1 2 2 1 2 2 2 1
## [186] 2 2 1 2 1 2 1 2 2 2 2 2 2 1 2 1 2 2 1 2 1 2 2 1 2 1 2 1 2 2 2 1 2 1 2 1 2
## [223] 1 2 2 2 2 2 2 2 2 1 2 1 2 1 1 2 2 1 2 1 2 1 2 2 1 2 1 2 1 2 2 2 2 2 2 2 1
## [260] 2 2 2 1 2 1 1 2 2 1 2 1 2
unique(dp4$clusterLabels) # Cantidad de clústeres para el conjunto de datos Faithful
## [1] 2 1

Dichas etiquetas de clúster son las resultantes a partir de la muestra final del algoritmo MCMC, según lo indican los autores del paquete. De esta forma, es posible graficar lo obtenido con anterioridad:

plot(dp4) # Gráfico del análisis multivariado de clústeres

Es importante mencionar que, para este caso en particular, resultaron únicamente dos clústeres de clasificación, sin embargo esta funcionalidad está generalizada para las dimensiones que sean necesarias (no hay límite en el número de agrupaciones).