R Code 2, Beta Binomial


#plot beta densities

a=1
b=1
p=seq(0,1,.002)
plot(p,dbeta(p,a,b),typ='l')

#Brute force approach to the binomial

y~Binomal(p,N)
p~beta(a,b)

find p|y

y=7
n=10
a=1
b=1

product=function(p) dbinom(y,n,p)*dbeta(p,a,b)
p.data=integrate(product,0,1)$value

plot(p,product(p)/p.data,typ="l",lwd=3,col="green",lty=2)
lines(p,dbeta(p,y+a,n-y+b),col='darkred')
lines(p,dbeta(p,a,b),col='darkblue')

Research Methods 3020