nma — Oct 16, 2013, 12:11 PM
function( DistPriori , vecDatos , credMass=0.95 , saveGr=FALSE ) {
# Inferencia Bayesiana para la verosimilitud de Bernoulli y distribucion a priori Beta.
# Argumentos de Entrada:
# DistPriori
# vector con los valores del parametro de la distribucion a priori beta.
# vecDatos
# la muestra es un vector de 1's y 0's.
# credMass
# probabilidad del IDMA.
# Salida:
# distPost
# vector con los valores del parametro de la distribucion a posteriori beta.
# Graficos:
# Crea un grafico de tres paneles con la a priori, verosimilitud
# y posteriori con el IDMA
# Ejemplo de uso:
# > distPost = BernBeta( DistPriori=c(1,1) , vecDatos=c(1,0,0,1,1) )
# Necesitas que R compruebe la funcion utilizar lo cual
# se hace mediante source("BernBeta.R"), tras nombrar su carpeta contenedora
# como directorio de trabajo getwd(trayectoria)
# Comprobacion de errores en los argumentos de entrada:
if ( length(DistPriori) != 2 ) {
stop("DistPriori debe tener dos componentes.") }
if ( any( DistPriori <= 0 ) ) {
stop("Las componentes de DistPriori deben ser positivas.") }
if ( any( vecDatos != 1 & vecDatos != 0 ) ) {
stop("vecDatos debe ser un vector de 1s y 0s.") }
if ( credMass <= 0 | credMass >= 1.0 ) {
stop("credMass debe estar entre 0 y 1.") }
# Renombrar los parametros de la distribucion a priori, por conveniencia:
a = DistPriori[1]
b = DistPriori[2]
# Se crean ciertos resumenes de los datos:
z = sum( vecDatos == 1 ) # numbero de 1's en vecDatos
N = length( vecDatos ) # numbero de individuos en vecDatos
# Se calculan los parametros de la distribucion a posteriori:
distPost = c( a+z , b+N-z )
# Calcula, p(D):
pData = beta( z+a , N-z+b ) / beta( a , b )
# Determinar los limites del IDMA.
# Llamada a la funcion HDIofICDF (debemos tenerlo nosotros).
source( "HDIofICDF.R" )
hpdLim = HDIofICDF( qbeta , shape1=distPost[1] , shape2=distPost[2] , credMass=credMass )
# Realiza los graficos:
# Construye la rejilla de valores de theta, utiles para realizar graficos
amplCaja = 0.005 # Valor pequeño para posicionar Theta en grafico.
Theta = seq( from = amplCaja/2 , to = 1-(amplCaja/2) , by = amplCaja )
# Calcula la priori de cada valor de theta
pTheta = dbeta( Theta , a , b )
# Calcula la verosimilitud de los datos en cada valor de theta.
pDataGivenTheta = Theta^z * (1-Theta)^(N-z)
# Calcula la posteriori de los datos en cada valor de theta.
pThetaGivenData = dbeta( Theta , a+z , b+N-z )
# Abre una ventana con tres paneles.
source("openGraphSaveGraph.R") # read in graph functions
openGraph(width=7,height=10,mag=0.7) # open a window for the graph
layout( matrix( c( 1,2,3 ) ,nrow=3 ,ncol=1 ,byrow=FALSE ) ) # 3x1 panels
par( mar=c(3,3,1,0) , mgp=c(2,1,0) , mai=c(0.5,0.5,0.3,0.1) ) # margin specs
maxY = max( c(pTheta,pThetaGivenData) ) # max y for plotting
# Dibuja la priori.
plot( Theta , pTheta , type="l" , lwd=3 ,
xlim=c(0,1) , ylim=c(0,maxY) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote(p(theta)) , cex.lab=1.5 ,
main="Prior" , cex.main=1.5 , col="skyblue" )
if ( a > b ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.0*max(pThetaGivenData) ,
bquote( "beta(" * theta * "|" * .(a) * "," * .(b) * ")" ) ,
cex=2.0 ,adj=textadj )
# Dibuja la verosimilitud: p(data|theta)
plot( Theta , pDataGivenTheta , type="l" , lwd=3 ,
xlim=c(0,1) , cex.axis=1.2 , xlab=bquote(theta) ,
ylim=c(0,1.1*max(pDataGivenTheta)) ,
ylab=bquote( "p(D|" * theta * ")" ) ,
cex.lab=1.5 , main="Likelihood" , cex.main=1.5 , col="skyblue" )
if ( z > .5*N ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.0*max(pDataGivenTheta) , cex=2.0 ,
bquote( "Data: z=" * .(z) * ",N=" * .(N) ) ,adj=textadj )
# Dibuja la posteriori.
plot( Theta , pThetaGivenData ,type="l" , lwd=3 ,
xlim=c(0,1) , ylim=c(0,maxY) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote( "p(" * theta * "|D)" ) ,
cex.lab=1.5 , main="Posterior" , cex.main=1.5 , col="skyblue" )
if ( a+z > b+N-z ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.00*max(pThetaGivenData) , cex=2.0 ,
bquote( "beta(" * theta * "|" * .(a+z) * "," * .(b+N-z) * ")" ) ,
adj=textadj )
text( textx , 0.75*max(pThetaGivenData) , cex=2.0 ,
bquote( "p(D)=" * .(signif(pData,3)) ) , adj=textadj )
# Marca el IDMA en la posteriori.
hpdHt = mean( c( dbeta(hpdLim[1],a+z,b+N-z) , dbeta(hpdLim[2],a+z,b+N-z) ) )
lines( c(hpdLim[1],hpdLim[1]) , c(-0.5,hpdHt) , type="l" , lty=2 , lwd=1.5 )
lines( c(hpdLim[2],hpdLim[2]) , c(-0.5,hpdHt) , type="l" , lty=2 , lwd=1.5 )
lines( hpdLim , c(hpdHt,hpdHt) , type="l" , lwd=2 )
text( mean(hpdLim) , hpdHt , bquote( .(100*credMass) * "% IDMA" ) ,
adj=c(0.5,-1.0) , cex=2.0 )
text( hpdLim[1] , hpdHt , bquote(.(round(hpdLim[1],3))) ,
adj=c(1.1,-0.1) , cex=1.2 )
text( hpdLim[2] , hpdHt , bquote(.(round(hpdLim[2],3))) ,
adj=c(-0.1,-0.1) , cex=1.2 )
# Construye nombrearchivo para el grafico guardado, y lo guarda.
if ( saveGr ) {
filename = paste( "BernBeta_",a,"_",b,"_",z,"_",N ,sep="")
saveGraph( file = nombrearchivo , type="eps" )
}
return( distPost )
}
function( DistPriori , vecDatos , credMass=0.95 , saveGr=FALSE ) {
# Inferencia Bayesiana para la verosimilitud de Bernoulli y distribucion a priori Beta.
# Argumentos de Entrada:
# DistPriori
# vector con los valores del parametro de la distribucion a priori beta.
# vecDatos
# la muestra es un vector de 1's y 0's.
# credMass
# probabilidad del IDMA.
# Salida:
# distPost
# vector con los valores del parametro de la distribucion a posteriori beta.
# Graficos:
# Crea un grafico de tres paneles con la a priori, verosimilitud
# y posteriori con el IDMA
# Ejemplo de uso:
# > distPost = BernBeta( DistPriori=c(1,1) , vecDatos=c(1,0,0,1,1) )
# Necesitas que R compruebe la funcion utilizar lo cual
# se hace mediante source("BernBeta.R"), tras nombrar su carpeta contenedora
# como directorio de trabajo getwd(trayectoria)
# Comprobacion de errores en los argumentos de entrada:
if ( length(DistPriori) != 2 ) {
stop("DistPriori debe tener dos componentes.") }
if ( any( DistPriori <= 0 ) ) {
stop("Las componentes de DistPriori deben ser positivas.") }
if ( any( vecDatos != 1 & vecDatos != 0 ) ) {
stop("vecDatos debe ser un vector de 1s y 0s.") }
if ( credMass <= 0 | credMass >= 1.0 ) {
stop("credMass debe estar entre 0 y 1.") }
# Renombrar los parametros de la distribucion a priori, por conveniencia:
a = DistPriori[1]
b = DistPriori[2]
# Se crean ciertos resumenes de los datos:
z = sum( vecDatos == 1 ) # numbero de 1's en vecDatos
N = length( vecDatos ) # numbero de individuos en vecDatos
# Se calculan los parametros de la distribucion a posteriori:
distPost = c( a+z , b+N-z )
# Calcula, p(D):
pData = beta( z+a , N-z+b ) / beta( a , b )
# Determinar los limites del IDMA.
# Llamada a la funcion HDIofICDF (debemos tenerlo nosotros).
source( "HDIofICDF.R" )
hpdLim = HDIofICDF( qbeta , shape1=distPost[1] , shape2=distPost[2] , credMass=credMass )
# Realiza los graficos:
# Construye la rejilla de valores de theta, utiles para realizar graficos
amplCaja = 0.005 # Valor pequeño para posicionar Theta en grafico.
Theta = seq( from = amplCaja/2 , to = 1-(amplCaja/2) , by = amplCaja )
# Calcula la priori de cada valor de theta
pTheta = dbeta( Theta , a , b )
# Calcula la verosimilitud de los datos en cada valor de theta.
pDataGivenTheta = Theta^z * (1-Theta)^(N-z)
# Calcula la posteriori de los datos en cada valor de theta.
pThetaGivenData = dbeta( Theta , a+z , b+N-z )
# Abre una ventana con tres paneles.
source("openGraphSaveGraph.R") # read in graph functions
openGraph(width=7,height=10,mag=0.7) # open a window for the graph
layout( matrix( c( 1,2,3 ) ,nrow=3 ,ncol=1 ,byrow=FALSE ) ) # 3x1 panels
par( mar=c(3,3,1,0) , mgp=c(2,1,0) , mai=c(0.5,0.5,0.3,0.1) ) # margin specs
maxY = max( c(pTheta,pThetaGivenData) ) # max y for plotting
# Dibuja la priori.
plot( Theta , pTheta , type="l" , lwd=3 ,
xlim=c(0,1) , ylim=c(0,maxY) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote(p(theta)) , cex.lab=1.5 ,
main="Prior" , cex.main=1.5 , col="skyblue" )
if ( a > b ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.0*max(pThetaGivenData) ,
bquote( "beta(" * theta * "|" * .(a) * "," * .(b) * ")" ) ,
cex=2.0 ,adj=textadj )
# Dibuja la verosimilitud: p(data|theta)
plot( Theta , pDataGivenTheta , type="l" , lwd=3 ,
xlim=c(0,1) , cex.axis=1.2 , xlab=bquote(theta) ,
ylim=c(0,1.1*max(pDataGivenTheta)) ,
ylab=bquote( "p(D|" * theta * ")" ) ,
cex.lab=1.5 , main="Likelihood" , cex.main=1.5 , col="skyblue" )
if ( z > .5*N ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.0*max(pDataGivenTheta) , cex=2.0 ,
bquote( "Data: z=" * .(z) * ",N=" * .(N) ) ,adj=textadj )
# Dibuja la posteriori.
plot( Theta , pThetaGivenData ,type="l" , lwd=3 ,
xlim=c(0,1) , ylim=c(0,maxY) , cex.axis=1.2 ,
xlab=bquote(theta) , ylab=bquote( "p(" * theta * "|D)" ) ,
cex.lab=1.5 , main="Posterior" , cex.main=1.5 , col="skyblue" )
if ( a+z > b+N-z ) { textx = 0 ; textadj = c(0,1) }
else { textx = 1 ; textadj = c(1,1) }
text( textx , 1.00*max(pThetaGivenData) , cex=2.0 ,
bquote( "beta(" * theta * "|" * .(a+z) * "," * .(b+N-z) * ")" ) ,
adj=textadj )
text( textx , 0.75*max(pThetaGivenData) , cex=2.0 ,
bquote( "p(D)=" * .(signif(pData,3)) ) , adj=textadj )
# Marca el IDMA en la posteriori.
hpdHt = mean( c( dbeta(hpdLim[1],a+z,b+N-z) , dbeta(hpdLim[2],a+z,b+N-z) ) )
lines( c(hpdLim[1],hpdLim[1]) , c(-0.5,hpdHt) , type="l" , lty=2 , lwd=1.5 )
lines( c(hpdLim[2],hpdLim[2]) , c(-0.5,hpdHt) , type="l" , lty=2 , lwd=1.5 )
lines( hpdLim , c(hpdHt,hpdHt) , type="l" , lwd=2 )
text( mean(hpdLim) , hpdHt , bquote( .(100*credMass) * "% IDMA" ) ,
adj=c(0.5,-1.0) , cex=2.0 )
text( hpdLim[1] , hpdHt , bquote(.(round(hpdLim[1],3))) ,
adj=c(1.1,-0.1) , cex=1.2 )
text( hpdLim[2] , hpdHt , bquote(.(round(hpdLim[2],3))) ,
adj=c(-0.1,-0.1) , cex=1.2 )
# Construye nombrearchivo para el grafico guardado, y lo guarda.
if ( saveGr ) {
filename = paste( "BernBeta_",a,"_",b,"_",z,"_",N ,sep="")
saveGraph( file = nombrearchivo , type="eps" )
}
return( distPost )
}