domingo, 23 de octubre de 2011

Si es Bayes, es bueno...


Los métodos de inferencia vistos hasta el momento forman parte de lo que se da en llamar “inferencia frecuentista”. El modelo subyacente es que, verificando ciertos supuestos, es posible ajustar una distribución (normal u otra) a mis datos para inferir conclusiones acerca de la población.
La inferencia bayesiana tiene otra visión. En general parte de que la distribución “real” de la población es algo imposible de conocer. En lugar de tomar una distribución como la verdadera, se toma lo que conoce “a priori” el investigador para elegir una distribución que se ajuste a ello y trabajar inicialmente bajo ese supuesto. Luego se toma una muestra de la población y con los datos obtenidos de ella se realiza un “ajuste” de la distribución empírica, llevando a una distribución que podría ser diferente a la anterior. Es lo que se conoce como distribución “a posteriori”.
O sea que, proponemos un modelo de acuerdo a lo que conocemos, tomamos evidencia para contrastarla contra el modelo y corregimos el modelo de acuerdo a lo aprendido.

Para tener un ejemplo intentaré estimar cuál es la proporción de estudiantes de una universidad que duerme al menos 8 horas.
Este ejemplo está tomado por completo del libro “Bayesian computation with R” de Jim Albert.

Supongamos que un investigador está interesado en conocer esa proporción. El verdadero valor de la proporción le es desconocido.
De acuerdo a la información de otras publicaciones y experiencia previa, piensa 2 cosas:
·         La mejor predicción que puede hacer es de 30%
·         Es muy probable que se halle por debajo del 50%

Construimos un rango de valores posibles para la proporción. Llamaremos “p” a la proporción. En R escribimos:

> p=seq(0.05,0.95,by=0.1)
> p
 [1] 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95

Esta instrucción construye un rango entre 0.05 (5%) y 0.95 (95%) que aumenta de a 0.1 (10%). Es un rango de valores posibles de p que nos servirá para construir una distribución a priori.
Ahora a cada uno de estos valores le otorgamos un puntaje entre 0 y 10 de acuerdo a lo que pensamos (esta escala es arbitraria). Estos puntajes son llamados “pesos”.

> pesos = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
> pesos
[1] 1.0 5.2 8.0 7.2 4.6 2.1 0.7 0.1 0.0 0.0

Y convertimos esos pesos a un formato de probabilidad. Para hacerlo, dividimos cada peso por la suma total de los mismos. A estas probabilidades se las llama “priors” o probabilidades “a priori”.

> prior = pesos/sum(pesos)
> prior

Esto nos permite crear una distribución que indica una probabilidad para cada valor de p. Para visualizarla podemos usar un gráfico:

> plot(p, prior, type = "h")

El gráfico obtenido nos muestra nuestra distribución “a priori”. Ahora la tarea será obtener una muestra de estudiantes. Lo hacemos y obtenemos datos de 27 estudiantes. Luego de revisar las encuestas notamos que hay 11 estudiantes que duermen más de 8 horas. Esto nos da un número de “éxitos” (11) y “fracasos” (16). Colocamos esos valores dentro de un objeto:

> data=c(11,16)
> data
[1] 11 16

El uso de la instrucción “c” es para concatenar varios valores y en este caso, asignarlos a un solo objeto.
Ahora hay que cargar un paquete de R llamado “LearnBayes”. Por supuesto, tiene que estar previamente instalado.

> library("LearnBayes")

Si no estuviese instalado nos daría un error diciendo que no encuentra el paquete. Procedemos como en el post que explica sobre instalación de paquetes de R y lo instalamos. Luego volvemos a tipear la instrucción anterior.

Luego de cargado el paquete escribimos el siguiente comando:

> post = pdisc(p, prior, data)
> post

Esto nos da las probabilidades a posteriori. Es bastante claro que son diferentes a las probabilidades a priori.
Para visualizarlas en forma comparativa escribimos una serie de instrucciones que no explicaré en detalle, pero nos permiten ver una gráfica sobre la otra. En caso de que a alguien le interese la explicación detallada me lo puede pedir. No obstante, se puede obtener ayuda de R escribiendo “?” seguido de la instrucción, o help(“instrucción”).

> library(lattice)
> PRIOR=data.frame("prior",p,prior)
> POST=data.frame("posterior",p,post)
> names(PRIOR)=c("Type","P","Probability")
> names(POST)=c("Type","P","Probability")
> data=rbind(PRIOR,POST)
>xyplot(Probability~P|Type,data=data,layout=c(1,2),type="h",lwd=3,col="black")

Esto nos muestra que las 2 distribuciones, la a priori construida con las suposiciones del investigador y la a posteriori, calculada usando la a priori y corrigiendo de acuerdo a las observaciones, son diferentes.

Existen funciones matemáticas para calcular las distribuciones a posteriori a partir de las a priori. Estas funciones son diferentes en cada situación y mi objetivo en este post no es ser extenso en éste tema.
Baste decir, que con esta metodología es posible estimar distribuciones pero también parámetros, como medias, sd y otros.

Con estos conocimientos podemos ir, ahora sí, a continuar con el aprendizaje del paquete PhViD para detección de señales.



No hay comentarios:

Publicar un comentario en la entrada