Created
September 3, 2025 06:53
-
-
Save abikoushi/5a098242ce0eb8f8a90f0b2d1d350096 to your computer and use it in GitHub Desktop.
Highest density interval of beta distribution
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # the following are used as reference | |
| # https://stats.stackexchange.com/questions/381520/how-can-i-estimate-the-highest-posterior-density-interval-from-a-set-of-x-y-valu | |
| hdi_beta = function(shape1, shape2, coverage, n){ | |
| x = seq(0,1, length.out=n) | |
| best = 0.0 | |
| for (ai in 1 : (length(x) - 1)){ | |
| for (bi in (ai + 1) : length(x)){ | |
| mass = pbeta(x[bi], shape1, shape2) - pbeta(x[ai], shape1, shape2) | |
| if (mass >= coverage && (mass / (x[bi] - x[ai])) > best){ | |
| best = mass / (x[bi] - x[ai]) | |
| ai.best = ai | |
| bi.best = bi | |
| } | |
| } | |
| } | |
| c(x[ai.best], x[bi.best]) | |
| } | |
| res = hdi_beta(3, 2, 0.95, 1001) | |
| diff(pbeta(res,3,2)) | |
| curve(dbeta(x,3,2)) | |
| abline(v = res, lty = 2) | |
| abline(h = dbeta(res,3,2), lty=2) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment