###--------------------------------------------------------------------------------------### ### Éder ###--------------------------------------------------------------------------------------### ### Solução analitica, númerica e por simulação do modelo # X ~ B(n,p) # p ~ Beta(alfa,beta) ###------------------------------------------------------------### ###------------------------------------------------------------### require(sfsmisc) require(latticeExtra) require(MASS) #browseURL('http://cs.illinois.edu/class/sp10/cs598jhm/Slides/Lecture02HO.pdf') ###------------------------------------------------------------### ###------------------------------------------------------------### ### grid de p p <- seq(0,0.99999,by=0.001) ### Priori alfa <- 10 beta <- 10 p.priori <- dbeta(p,alfa,beta) ### Verossimilhança n <- 1000 x <- rbinom(1,n,0.3) vero <- function(p,n,x){exp(sum(dbinom(x,n,p,log=TRUE)))} p.vero <- apply(matrix(p),1,vero,n=n,x=x) ###------------------------------------------------------------### ###------------------------------------------------------------### ### Solução analitica ### Posteriori p.posteA <- dbeta(p,alfa+sum(x),beta+sum(n-x)) ### Plotando doubleYScale(xyplot(p.priori + p.posteA ~ p, type = "l",lwd=3), xyplot(p.vero ~ p, type = "l",lwd=2,lty=2), style1 = 0, style2 = 3, add.ylab2 = TRUE, text = c("Priori", "Posteriori", "Verossimilhança"), columns = 3) ### confirmando se a posteriori é uma fdp integrate.xy(p,p.posteA) ###------------------------------------------------------------### ###------------------------------------------------------------### ### Integração númerica para normalização ### posteriori p.posteN <- (p.priori*p.vero)/(integrate.xy(p,p.priori*p.vero)) ###Quadratura Gaussiana N <- 100 nos <- ((1-(0))*(1:N-0.5))/N nos Ia <- (1-(0)) * sum(apply(as.matrix(nos),1,function(p){vero(p,n,x)*dbeta(p,alfa,beta)}))/N Ia (integrate.xy(p,p.priori*p.vero)) ### Plotando doubleYScale(xyplot(p.priori + p.posteN ~ p, type = "l",lwd=2), xyplot(p.vero ~ p,type = "l",lwd=2,lty=2), style1 = 0, style2 = 3, add.ylab2 = TRUE, text = c("Priori", "Posteriori", "Verossimilhança"), columns = 3) ### confirmando se a posteriori é uma fdp integrate.xy(p,p.posteN) ###------------------------------------------------------------### ###------------------------------------------------------------### ### Amostragem da posteriori ns <- 10000 theta_chapeu <- sum(x)/(n*length(x)) theta_i <- rbeta(ns,alfa,beta) u_i <- runif(ns,0,1) crite <- u_i <= ((dbeta(theta_i,alfa,beta)*apply(matrix(theta_i),1,vero,n=n,x=x))/ (dbeta(theta_chapeu,alfa,beta)*vero(theta_chapeu,n=n,x=x))) a.posteriori <- theta_i[crite] mean(a.posteriori,na.rm=TRUE) ### Taxa Aceitação sum(crite)/ns ###------------------------------------------------------------### ###------------------------------------------------------------### ### Comparando os resultados hist(a.posteriori,prob=TRUE) rug(a.posteriori) lines(density(a.posteriori)) lines(p,p.posteA,col='red',lwd=3) lines(p,p.posteN,col='blue',lty=2) legend('topleft',c('Amostragem','Analitico','Númerica'),lty=c(1,1,2),col=c('black','red','blue')) ### Intervalos via verosimilhança aproximado theta_chapeu+c(-1,1)*1.96*sqrt((theta_chapeu*(1-theta_chapeu))/n) ### IC amostragem quantile(a.posteriori,c(0.025,0.975)) ### Analitico da conjugada qbeta(c(0.025,0.975),alfa+sum(x),beta+sum(n-x)) ###--------------------------------------------------------------------------------------### ### FIM... ###--------------------------------------------------------------------------------------###