累積密度関数の逆関数を与えてHighest Density Intervalを算出する

HDI=\(InvCDF=qbeta,cred=0.95,...){
opt=optimize(\(p) InvCDF(p+cred,...) - InvCDF(p,...),c(0,1-cred))
lwr=InvCDF(opt$min,...)
upr=lwr+opt$obj
c(lwr,upr)
}

ガンマ分布で検証
lu=HDI(qgamma,shape=20,rate=24)
lu
curve(dgamma(x,shape=20,rate=24),0,2)
pgamma(lu[2],shape=20,rate=24)-pgamma(lu[1],shape=20,rate=24)
diff(lu) < qgamma(0.975,20,24)-qgamma(0.025,20,24)


ベータ分布で検証
lu=HDI(qbeta,shape1=20,shape2=24)
lu
curve(dbeta(x,20,24),0,1)
pbeta(lu[2],20,24)-pbeta(lu[1],20,24)
diff(lu) < qbeta(0.975,20,24)-qbeta(0.025,20,24)