Pensamientos y Libros

22 diciembre 2006

Instalando JAGS para estimación Bayesiana.

Filed under: Estadística Bayesiana,Matemáticas — Toni Cebrián @ 12:03

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()

Representación de el parámetro phi obtenido en las simulaciones.

Dejar un comentario »

Aún no hay comentarios.

RSS feed for comments on this post. TrackBack URI

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

Blog de WordPress.com.

A %d blogueros les gusta esto: