path <- commandArgs()[6] ############################## pb <- function(p, a, b) { exp((a-1) * log(p) + (b-1) * log(1-p) - lbeta(a, b)) } pv <- seq(0, 1, by = 0.0025) z1 <- pb(pv, 1, 2/3)[!is.na(pb(pv, 1, 2/3)) & !is.infinite(pb(pv, 1, 2/3))] p1 <- pv[!is.na(pb(pv, 1, 2/3)) & !is.infinite(pb(pv, 1, 2/3))] z2 <- pb(pv, 1 + 650, 2/3 + 350)[!is.na(pb(pv, 1 + 650, 2/3 + 350)) & !is.infinite(pb(pv, 1 + 650, 2/3 + 350))] p2 <- pv[!is.na(pb(pv, 1 + 650, 2/3 + 350)) & !is.infinite(pb(pv, 1 + 650, 2/3 + 350))] (1 + 650) / (1 + 2/3 + 1000) (1 + 650) * (2/3 + 350) / ((2/3 + 1 + 1000)^2 * (2/3 + 1 + 1000 + 1)) ######## prior ######## postscript( file=paste(path, "Ex09-1.eps", sep=""), width = 9, height = 9, pointsize = 24, horizontal = FALSE, onefile = FALSE, paper = "special", family="Palatino" ) par(mar=c(3.2, 3.2, 1.5, 0.5), mgp=c(2.0, 0.7, 0)) par(bty="l") plot( p1, z1, xlim = c(0, 1), type = 'l', lty = 1, col = "#E64B6B", lwd = 6, ylim = c(0, 5), xlab = expression(theta), ylab = expression(paste("p(", theta, ")")), main = "prior density", panel.first = grid(NA, NULL, lty = 2, col = "#E9DECA") ) text( 0, max(z1) * 0.95, expression(paste("Beta(1, 2/3)")), pos = 4, offset = 0, cex=.8 ) dev.off() ######## posterior ######## postscript( file=paste(path, "Ex09-2.eps", sep=""), width = 9, height = 9, pointsize = 24, horizontal = FALSE, onefile = FALSE, paper = "special", family="Palatino" ) par(mar=c(3.2, 3.2, 1.5, 0.5), mgp=c(2.0, 0.7, 0)) par(bty="l") plot( p2, z2, xlim = c(0, 1), type = 'l', lty = 1, col = "#E64B6B", lwd = 6, xlab = expression(theta), ylab = expression(paste("p(", theta, " |", y, ")")), main = "posterior density", panel.first = grid(NA, NULL, lty = 2, col = "#E9DECA") ) text( 0, max(z2) * 0.95, expression(paste("Beta(1 + 650, 2/3 + 350)")), pos = 4, offset = 0, cex=.8 ) dev.off()