placenta.r_function (pr.inf=c(0,0,10,100,1000,10000)) 
{
# Bayesian estimation of a probability from binomial data
# Example at pag. 39, sec 2.5, Gelman reference text

# Analysis using a uniform prior [Beta(alpha=1,beta=1)]
# as well as several Beta priors which are increasingly
# concentrated around the mean set at 0.485, the known
# mean of the general population (female births)

# After fixed different values of alpha+beta (whence values
# of alpha, beta result accordingly with the assigned prior 
# mean), the code draws posteriors, likelihoods and priors 
# for each case.

# To evaluate the sensitivity of posterior inference to 
# the proposed prior, the 95% posterior intervals of theta
# and the general population mean (0.485) are drawn.

# Finally, a summary table of the analysis is printed to a 
# a tex file "placenta.tex".
  
alpha.beta_pr.inf+2
alpha_c(1,.485*alpha.beta[-1])
beta_c(1,alpha.beta[-1]-alpha[-1])

layout(matrix(1:6,3,2,byrow=T))
par(mar=c(2,2,2,2))

theta_seq(.35,.55,.001)
posum_c()

for(i in 1:length(alpha))
 {
  plot(theta,dbeta(theta,437+alpha[i],543+beta[i]),
  type="n",bty="l")
  text(x=.4,
  y=dbeta((437+alpha[i])/(980+alpha[i]+beta[i]),
  437+alpha[i],543+beta[i])-2.5,
  paste("a+b-2=",as.character(pr.inf[i])),cex=1.5)

# likelihood. It has been scaled (per 1000) to be 
# comparable with posterior at alpha+beta=2.

  lines(theta,1000*choose(980,437)*theta^437*(1-theta)^543,
  lwd=2,col="red")

# prior

  lines(theta,dbeta(theta,alpha[i],beta[i]),
  lwd=2,col="aquamarine")

# posterior

  lines(theta,dbeta(theta,alpha[i]+437,beta[i]+543),
  lwd=2,col="orange")

# posterior inference

  ci_qbeta(c(.025,.975),437+alpha[i],543+beta[i])
  abline(v=c(0.485))
  abline(v=ci,col="orange")
  posum_rbind(posum,c(pr.inf[i],alpha[i]/(alpha[i]+beta[i]),
  round(qbeta(.5,alpha[i]+437,beta[i]+543),3),
  paste("[",as.character(round(ci[1],3)),",",
  as.character(round(ci[2],3)),"]")) )
  }
tab_xtable(posum, align=c("r","r","c","c","c"),
digits=c(0,0,3,3,3),caption="A former sensitivity analysis")
print.xtable(tab,file="placenta.tex")
return(posum)
}