Python中的分析最高密度区间（最好用于Beta发行版）

Question

我想知道是否有人知道可靠和快速的分析HDI计算，最好是针对贝塔函数

人类发展指数的定义被称为“最高后密度区”

我正在寻找具有以下I/O的函数：

输入

• credMass-可信区间质量（例如，95%可信区间的{}）
• a-形状参数（例如，硬币投掷的头数）
• b-形状参数（例如，抛硬币的尾巴数）

输出

• ci_min-质量的最小界限（介于{}到{}之间的值）
• ci_max-质量的最大界限（介于{}到{}之间的值）

一种方法是加速我发现的这个脚本 在the same said question（从R改编为Python）中，摘自John K.Kruschke的《做贝叶斯数据分析》一书）。我使用过这个解决方案，它非常可靠，但是对于多个调用来说它有点太慢了。加速100倍甚至10倍会非常有帮助

from scipy.optimize import fmin
from scipy.stats import *

def HDIofICDF(dist_name, credMass=0.95, **args):
    # freeze distribution with given arguments
    distri = dist_name(**args)
    # initial guess for HDIlowTailPr
    incredMass =  1.0 - credMass

    def intervalWidth(lowTailPr):
        return distri.ppf(credMass + lowTailPr) - distri.ppf(lowTailPr)

    # find lowTailPr that minimizes intervalWidth
    HDIlowTailPr = fmin(intervalWidth, incredMass, ftol=1e-8, disp=False)[0]
    # return interval as array([low, high])
    return distri.ppf([HDIlowTailPr, credMass + HDIlowTailPr])

用法

print HDIofICDF(beta, credMass=0.95, a=5, b=4)

警告！一些解决方案将此HDI与等尾间隔解决方案（该问题称之为“中央可信区域”）混淆，后者更容易计算，但不能回答相同的问题。（例如，见Kruschke的Why HDI and not equal-tailed interval?）

此外，这个问题并不涉及我在PyMC3中看到的MCMC方法（pymc3.stats.hpd(a)，其中a是一个随机变量样本），而是涉及一个分析解

谢谢大家!