Tabla de contenidos
- Introducción y motivación
- Instalación
- Comprobando que todo está en orden y funciona
- Ejemplo de uso
Introducción y motivación
Ultimamente he estado leyendo los trabajos de Radford M. Neal y de David J.C. MacKay acerca de las Redes Neuronales Bayesianas y de los Procesos Gaussianos y me ha picado la curiosidad sobre la estimación bayesiana. También han servido para avivar el gusanillo bayesiano el libro Scientific Reasoning: The Bayesian Approach que leí hace poco y el que pretendo leerme en breve Probability Theory: The Logic of Science.
Con todo me puse a investigar cómo podía ponerme a trastear con esto de la inferencia bayesiana en Linux y encontré el programita JAGS para inferencia bayesiana utilizando el método Monte Carlo utilizando Cadenas de Markov (MCCM). Con eso y con la amable ayuda de su autor, Martyn Plummer, he conseguido tener funcionando el invento. Espero que este sea el primer artículo donde vaya contando mis experiencias al respecto.
Instalación
Lo primero que tenemos que hacer es instalar en nuestra máquina Linux las librerías de desarrollo y los paquetes de los que depende JAGS. Para ello instalamos los siguientes paquetes:
$> su
Password:
$> urpmi gcc-gfortran
$> urpmi libblas1.1-devel
$> urpmi liblapack3-devel
descargarnos el software de su página web, JAGS. Nos descargamos la última versión que exista, en nuestro caso la 0.97.1. De nuevo como usuario sin privilegios de root descomprimimos el software, configuramos las dependencias, compilamos e instalamos JAGS.
$> tar xvfz JAGS-0.97.tar.gz
$> ./configure
$> make
$> su
Password:
$> make install
Comprobando que todo está en orden y funciona
Ya tenemos instalado JAGS y nos queremos poner manos a la obra. Aunque JAGS puede funcionar de forma autónoma nos será mucho más fácil utilizarlo dentro del paquete estadístico R así que nos lo instalamos también:
$> urpmi R-base
Para poder utilizar el paquete rjags es preciso que previamente estén instalados en el entorno de R los paquetes coda y lattice. La forma más cómoda es que estos paquetes se instalen como root de forma que estén disponibles para el resto de usuarios (que me perdonen los seguidores de la seguridad paranoica). Desde R aquellos paquetes que estén en algún repositorio de CRAN se instalan de forma muy fácil. Primero se selecciona el espejo de CRAN desde donde descargar los paquetes y después con el comando de instalar él sólo hará la instalación:
$> R
> chooseCRANmirror()
Selección: 45
> install.packages("coda")
En el cuadro anterior se han omitido los mensajes intermedios.
En este momento simplemente escribiendo el siguiente comando en el entorno R y comprobando que no se produce ningún mensaje de error tenemos el entorno configurado y listo para usar.
> library(rjags)
Ejemplo de uso
Para ilustrar su funcionamiento realizaremos un pequeño ejemplo. Mediante Perl implementamos un pequeño generador de números aleatorios binarios. El comando utilizado puede verse más abajo.
$> perl -e 'for($i=0;$i<1000;$i+=1){print int(rand()+.5);print "n"}' | grep 1 | wc -l
que nos genera un listado de unos y ceros generados teóricamente por una distribución de Bernoulli de parámetro p=0.5 que después utilizaremos para obtener el número de 1s obtenidos. Lo que pretendemos hacer en este ejemplo práctico es tratar de estimar p desde un punto de vista Bayesiano y poder sacar conclusiones sobre la calidad de la generación de números aleatorios con el procedimiento anterior.
Ejecutando la línea de Perl me sale un valor de 529 para el experimento (cada ejecución dará un valor distinto). Por tanto sólo nos queda crear el modelo con la syntaxis BUGS,
model
{
phi ~ dunif(0,1);
y ~ dbin(phi,N);
}
y crear un archivo con los datos necesarios para el problema.
datos <- list(y=529,N=1000)
Después únicamente nos queda arrancar R, cargar la libreria rjags y empezar a interrogar al modelo.
>library(rjags)
Loading required package: coda
Loading required package: lattice
loading JAGS modules
basefunctions
baserngs
basesamplers
bugs
> source("binario.dat")
> modelo <- jags.model("binario.bugs",datos)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Checking graph
Graph Size: 5> modelo$update(1000)
Updating 1000
---------------------------------------| 1000
**************************************** 100%
> x <- model.samples(modelo, variable.names=c("phi"),n.iter=1000)
Updating 1000
---------------------------------------| 1000
**************************************** 100%
> summary(x$phi)
V1
Min. :0.4847
1st Qu.:0.5172
Median :0.5284
Mean :0.5285
3rd Qu.:0.5400
Max. :0.5812
donde puede verse que el valor más probable está en torno a 0.52, para el parámetro p.
Sin embargo, en nuestra estimación hemos utilizado una distribución uniforme en el rango de posibles valores de p. Esto no es una hipótesis muy realista pues tenemos que pensar que los programadores que hicieron la función rand() de Perl, tienen los conocimientos de programación y estadística suficientes como para esperar que en la implementación que hicieran el valor de p no difiera mucho de 0.5. Podemos probar a cambiar la distribución inicial de p para que su densidad inicial sea una distribución Normal concentrada en torno a 0.5 y por ejemplo su desviación estándar, sigma, valga 0.05, lo que nos asegura que el 68% de la masa de la función de densidad estará en el intervalo [0.45,0,55].
Creamos el nuevo modelo especificando una distribución normal y teniendo en cuenta que el parámetro tau de la distribución normal se llama precisión y viene dado por 1/sqrt(sigma), donde sqrt es la función raíz cuadrada. Por tanto el nuevo modelo será:
model{
phi ~ dnorm(0.5,400) T(0,1);
y ~ dbin(phi,N);
}
Es importante para los cálculos de JAGS que se especifique que el valor de phi debe truncarse a un valor entre 0 y 1 para que pueda considerarse un valor legal como parámetro de la distribución binomial. Esto lo hemos hecho añadiendo detrás de la distribución el parámetro de truncamiento, T(0,1).
Como antes, seguimos el proceso de compilar, actualizar e interrogar al model y como paso final obtenemos una representación gráfica de los datos obtenidos como estimación del parámetro phi.
> m <- jags.model("binario-normal.bugs",datos)Compiling model graph
Resolving undeclared variables
Allocating nodes
Checking graph
Graph Size: 7> m$update(1000)
Updating 1000
---------------------------------------| 1000
**************************************** 100%
> x <- model.samples(m,variable.name=c("phi"),n.iter=1000)
Updating 1000
---------------------------------------| 1000
**************************************** 100%
> x <- model.samples(m,variable.name=c("phi"),n.iter=1000)
> xfig()
> densityplot(x$phi)
> dev.close()
Pensamientos y Libros 