Аналитический интервал максимальной плотности в Python (предпочтительно для бета-версий)

Мне было интересно, знает ли кто-нибудь о надежном и быстром аналитическом расчете HDI, предпочтительно для бета-функций.

Определение HDI содержится в этом вопросе под названием Highest Posterior Область плотности.

Я ищу функцию, которая имеет следующий ввод-вывод:

Ввод

• credMass - масса достоверного интервала (например, 0.95 для 95% достоверного интервала)
• a - параметр формы (например, количество подбрасываний монеты HEADS)
• b - параметр формы (например, количество подбрасываний монеты TAILS)

вывод

• ci_min - минимальная граница массы (значение от 0. до ci_max)
• ci_max - максимальная граница массы (значение от ci_min до 1.)

Один из способов сделать это — ускорить этот скрипт, который я нашел в тот же самый вопрос (адаптированный из R для Python) и взятый из книги Джона К. Крушке «Выполнение байесовского анализа данных»). Я использовал это решение, которое достаточно надежно, но слишком медленно для нескольких вызовов. Ускорение в 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 с решением Equal Tail Interval (которое в указанном вопросе называется Центральным достоверным регионом), которое гораздо проще вычислить, но не дает ответа на тот же вопрос. (Например, см. статью Крушке Почему HDI, а не интервал с равными хвостами?< /эм>)