Visualización de la actualización bayesiana (y unas cuantas funciones de R)

2011-9-12 (Última modificación: 2011-9-12)

Me ha llegado noticia de una entrada en un blog, Visualizing Bayesian Updating, en el que se muestra visualmente cómo se actualiza la distribución a posteriori conforme aumenta el número de ensayos en un problema bayesiano simple. Explica también los fundamentos estadísticos del asunto.

Yo me limitaré a ofrecer una nueva versión del código —que no funcionaba copiando y pegando sin más— en el que he introducido ciertas modificaciones. Es el siguiente:

sim.bayes <- function(p=0.5, N=10, y.lim=15)
{
  plot( 1, xlim = c(0,1), ylim = c(0, y.lim),
        xlab = 'p', ylab = 'Posterior Density',
        type = "n" )
  legend('topright',
      legend=c('Prior','Updated Posteriors','Final Posterior'),
      lty=c(2,1,1),col=c('black','black','red'))

  exitos <- cumsum( rbinom( N, 1, p ) )

  foo <- function( exitos, n.pruebas, col = "black", lwd = 1, lty = 1 ){
    curve( dbeta( x, exitos + 1, n.pruebas - exitos + 1 ),
            add = TRUE, col = col, lwd = lwd, lty = lty )
    print( paste(exitos, "éxitos y ", n.pruebas - exitos, "fracasos") )
  }

  foo( 0, 0, lty = 2 )
  mapply( foo, exitos, 1:N )
  foo( exitos[N], N, col='red', lwd = 2 )
}

sim.bayes(p = 0.6, N = 20, y.lim = 5)

Quiero pensar que mis lectores encontrarán útil el ejemplo de uso de la función mapply (para recorrer dos vectores simultáneamente), curve (para representar gráficamente funciones).