1binomial.beta.mix=function(probs,betapar,data)
2{
3N=length(probs)
4s=data[1]; f=data[2]
5post.betapar=betapar+outer(rep(1,N),data)
6p=post.betapar[,1]/(post.betapar[,1]+post.betapar[,2])
7m.prob=exp(dbinom(s,size=s+f,prob=p,log=TRUE)+
8  dbeta(p,betapar[,1],betapar[,2],log=TRUE) -
9  dbeta(p,post.betapar[,1],post.betapar[,2],log=TRUE))
10
11post.probs=probs*m.prob/sum(probs*m.prob)
12return(list(probs=post.probs,betapar=post.betapar))
13}