nma — Sep 30, 2013, 12:50 PM
# Theta es el vector de valores posibles de theta (v.a.)
# nValsTheta es el número de valores posibles de theta (v.a.), es decir la dimensión del vector theta
nValsTheta = 3
# lo siguiente se hace para asignar los valores al vector Theta de forma homogenea
Theta = seq ( from = 1/(nValsTheta+1), to =nValsTheta/(nValsTheta+1), by =1/(nValsTheta+1))
# A priori:
pTheta = pmin(Theta,1-Theta)
pTheta = pTheta/sum(pTheta)
# Datos es un vector de correspondiente a una m.a.s. de Bernoulli's
Datos = c(1,1,1,0,0,0,0,0,0,0,0,0)
NCaras = sum(Datos == 1) #Binomial que indica número de caras en el total de intentos
NCruces = sum(Datos == 0)
# Verosimilitud:
pDatosCondicionadoTheta = (Theta^NCaras)*((1-Theta)^NCruces)
# Posteriori:
pDatos = sum(pDatosCondicionadoTheta*pTheta)
pThetaCondicionadoDatos = (pDatosCondicionadoTheta * pTheta) / pDatos # Teorema de Bayes
# Para dibujar los resultados
windows(7,10) #Crea una ventana con longitud y anchura determinada (en pulgadas)
layout(matrix(c(1,2,3), nrow=3, ncol=1, byrow=FALSE)) #Panel 3x1
par(mar=c(3,3,1,0)) #número de lineas para el margen: abajo, izquierda, arriba, derecha
par(mgp=c(2,1,0)) #número de lineas para leyendas
par(mai=c(0.5,0.5,0.3,0.1)) #longitud para el margen (pulgadas): abajo, izquierda, arriba, derecha
# Para dibujar la priori
plot(Theta, pTheta, type="h", lwd=3, main="Priori",
xlim=c(0,1), xlab=bquote(theta),
ylim=c(0,1.1*max(pThetaCondicionadoDatos)), ylab=bquote(p(theta)),
cex.axis=1.2, cex.lab=1.5, cex.main=1.5)
# Para dibujar la verosimilitud
plot(Theta, pDatosCondicionadoTheta, type="h", lwd=3, main="Verosimilitud",
xlim=c(0,1), xlab=bquote(theta),
ylim=c(0,1.1*max(pDatosCondicionadoTheta)), ylab=bquote(paste("p(D|",theta,")")),
cex.axis=1.2, cex.lab=1.5, cex.main=1.5)
text(.55,.85*max(pDatosCondicionadoTheta),cex=2.0,
bquote("D="*.(NCaras)*"C,"*.(NCruces)*"X"),adj=c(0,.5))
# Para dibujar la posteriori
plot(Theta, pThetaCondicionadoDatos, type="h", lwd=3, main="Posteriori",
xlim=c(0,1), xlab=bquote(theta),
ylim=c(0,1.1*max(pThetaCondicionadoDatos)), ylab=bquote(paste("p(",theta,"|D)")),
cex.axis=1.2, cex.lab=1.5, cex.main=1.5)
text(.55,.85*max(pThetaCondicionadoDatos),cex=2.0,
bquote("p(D)="*.(signif(pDatos,3))),adj=c(0,.5))