# Stat 241/541 fall 2014
# Example 3.4

# b(n,p,k) = (n choose k) p^k (1-p)^(n-k) = dbinom(k,size=n,prob=p)

# Enable simultaneous calculation for more than one prob:
BinProbs <- Vectorize(dbinom,vec="prob")


makePosterior <- function(p = c(0.45,0.5,0.55), n=10, prior = c(0.5,0.3,0.2)){

	bb <- BinProbs(0:n,size=n,prob=p)
	colnames(bb) <- p	# Now bb[i,j] = b(n,i,p_j)
	
	W <- bb %*% diag(prior)	# W[i,j] = b(n,i,p_j) x prior[j]
	Dprob <- apply(W,1,sum)		# Dprob[k] = prob(D_k)
	post <- diag(1/Dprob) %*% W	# makes each row sum to 1

# The last two matrix multiplications (using %*%) could be replaced 
#		applications of the sweep() function.

	Grid <- (0:n)/n
	matplot(Grid,post,
		type="l", 
		col=c("black","red","blue"),
		xlab="proportion of heads", 
		ylab = paste("posterior prob"),
		ylim=c(0,1),
		axes = F
	)	
	xvalues <- c(0,p,1)
	axis(1, at=xvalues, labels = as.character(xvalues)   )
	axis(2, at=c(0,1))

	title(main = paste("posterior from n = ", n, "tosses") )
	return(post)	# not really needed
}

draw.posterior = function(p = c(0.4,0.5,0.6),tosses = c(10,50,100),...){
	par(mfrow=c(length(tosses),1)) # stack the pictures
	for (n in tosses){
		makePosterior(p,n,...)
	}
	par(mfrow=c(1,1))
}