betabin.r_function (y,m,a,b) 
{
n_length(a)  # n=1 or >1&even
theta_seq(0.,1.0,.001)

upost_rbeta(10000,y+1,m-y+1)
uci_round(quantile(upost,prob=c(.025,.975),names=F),3)

maxupo_max(hist(upost,nclass=40,probability=T,plot=F)$density)
maxli_max(choose(m,y)*theta^y*(1-theta)^(m-y))

pr2po_matrix(,nrow=n,ncol=5)
# pr2poy_matrix(,nrow=n,ncol=4)

windows(canvas="black",gamma=2)
if(n>1)
 {layout(matrix(1:n,nrow=n/2,byrow=T))
  par(mar=c(2,2,2,2))}

for(i in 1:n)
{
if(a[i]==1 & b[i]==1) post_upost
else post_rbeta(10000,y+a[i],m-y+b[i])

ci_round(quantile(post,prob=c(.025,.975),names=F),3)
nci_round(ci.r(post),3)
logits_logit.r(post)
nlci_ci.r(logits)
nilci_round(c(invlogit.r(nlci[1]),invlogit.r(nlci[2])),3)

maxpo_max(hist(post,probability=T,plot=F)$density)
dmax_max(maxpo,maxupo)
hist(post,probability=T,nclass=40,col="yellow2",
     fg="white",col.axis="white",ylim=c(0,dmax))
lines(theta,maxupo/maxli*choose(m,y)*theta^y*(1-theta)^(m-y),
      lwd=2,col="red")
mtext(side=1,col="aquamarine",text=
  paste("a+b-2=",as.character(round(a[i]+b[i]-2,4)),
        "  a=",as.character(a[i])),cex=1)

abline(v=ci,col="magenta",lwd=2)
abline(v=nci,col="red",lwd=2)
abline(v=nilci,col="white",lwd=2)
abline(v=uci,col="green",lwd=2)

p_a[i]/(a[i]+b[i])
prv_round(p*(1-p)/(a[i]+b[i]+1),3)
pp_round((a[i]+y)/(a[i]+b[i]+m),3)
pov_round(pp*(1-pp)/(a[i]+b[i]+m+1),3)

pr2po[i,]_c(pp,ci,prv,pov)
}
return(pr2po)
}

# Internal calls in betabin.r are to the following functions:

logit.r_function (p) 
{
y_log(p/(1-p))
return(y)
}


invlogit.r_function (x) 
{
y_exp(x)/(1+exp(x))
return(y)
}

ci.r_function (champion) 
{
ci.mean_mean(champion)
ci.half_1.96*sd(champion)
return(c(ci.mean-ci.half,ci.mean+ci.half))
}