# Stat 241/541 fall 2005 # variations on the HHH vs. TTHH example setwd("/Users/pollard/Academic/Courses/241/241.fall2005/notes2005") posterior0 = function(p1=0.1, p2=0.5,p3=0.9,n=20){ k = 0:n L1 = p1^k * (1-p1)^(n-k) L2 = p2^k * (1-p2)^(n-k) L3 = p3^k * (1-p3)^(n-k) denom = L1 +L2 + L3 L1 = L1/denom L2 = L2/denom L3 = L3/denom plot(k,L1,pch=19,type="b") points(k,L2, pch="x",type="b") points(k,L3, pch=21,type="b") return(rbind(L1,L2,L3)) } posterior = function(p = c(0.1,0.5,0.9), n=20, prior = rep(1,n+1)){ kk = 0:n probs = outer(kk,p,function(k,p) {p^k *(1-p)^(n-k)}) probs = probs %*% diag(prior) totals = apply(probs, 1,sum) probs =sweep(probs, 1,totals,"/") matplot(kk/n,probs, type="b",pch=seq(p), xlab="k/n", ylab = paste("posterior prob given k heads"), axes = F) xvalues =c(0,p,1) axis(1, at=xvalues, labels = as.character(xvalues) ) axis(2, at=NULL) title(main = paste("posterior probabilities from ", n, "tosses"),tcl=1, ) return(probs) } draw.posterior = function(p = c(0.4,0.5,0.6),tosses = c(10,50,100),...){ old.par= par() par(mfrow=c(1,length(tosses))) for (n in tosses){ posterior(p,n,...) } par(mfrow=c(1,1)) }