beta HPD

October 17, 2013
By

(This article was first published on Xi'an's Og » R, and kindly contributed to R-bloggers)

While writing an introductory chapter on Bayesian analysis (in French), I came by the issue of computing an HPD region when the posterior distribution is a Beta B(α,β) distribution… There is no analytic solution and hence I resorted to numerical resolution (provided here for α=117.5, β=115.5):

f=function(p){

  # find the symmetric
  g=function(x){return(x-p*((1-p)/(1-x))^(115.5/117.5))}
  return(uniroot(g,c(.504,.99))$root)}

ff=function(alpha){

  # find the coverage
  g=function(x){return(x-p*((1-p)/(1-x))^(115.5/117.5))}
  return(uniroot(g,c(.011,.49))$root)}

and got the following return:

> ff(.95)
[1] 0.4504879
> f(ff(.95))
[1] 0.5580267

which was enough for my simple book illustration… Since (.450,558) is then the HPD region at credible level 0.95.


Filed under: Books, R, Statistics, Uncategorized, University life Tagged: beta distribution, book chapter, fREN, French paper, HPD region, pbeta(), R, uniroot()

To leave a comment for the author, please follow the link and comment on his blog: Xi'an's Og » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.