Блог

для стохастического моделирования мне нужно нарисовать множество случайных чисел, которые распределены бета-биномиально.

На данный момент я реализовал это таким образом (используя python):

 import scipy as scp
from scipy.stats import rv_discrete

class beta_binomial(rv_discrete):
       """
       creating betabinomial distribution by defining its pmf
       """
       def _pmf(self, k, a, b, n):
          return scp.special.binom(n,k)*scp.special.beta(k a,n-k b)/scp.special.beta(a,b)
  

таким образом, выборка случайного числа x может быть выполнена с помощью:

 betabinomial = beta_binomial(name="betabinomial")
x = betabinomial.rvs(0.5,0.5,3) # with some parameter 
  

Проблема в том, что выборка одного случайного числа занимает около 0,5 мс, что в моем случае доминирует над всей скоростью моделирования. Ограничивающим элементом является оценка бета-функций (или гамма-функций внутри них).

У кого-нибудь есть отличная идея, как ускорить выборку?

Что ж, вот рабочий и слегка протестированный код, который кажется более быстрым, используя свойство составного распределения бета-биномиального.

Мы делаем выборку p из бета-версии, а затем используем ее в качестве параметра для binomial. Если бы вы выбирали векторы большого размера, это было бы еще быстрее.

 import numpy as np

def sample_Beta_Binomial(a, b, n, size=None):
    p = np.random.beta(a, b, size=size)
    r = np.random.binomial(n, p)

    return r

np.random.seed(777777)
q = sample_Beta_Binomial(0.5, 0.5, 3, size=10)
print(q)
  

Результат равен

 [3 1 3 2 0 0 0 3 0 3]
  

Быстрый тест

 np.random.seed(777777)

n = 10
a = 2.
b = 2.
N = 100000

q = sample_Beta_Binomial(a, b, n, size=N)

h = np.zeros(n 1, dtype=np.float64) # histogram
for v in q: # fill it
    h[v]  = 1.0

h /= np.float64(N) # normalization
print(h)
  

выводит гистограмму

 [0.03752 0.07096 0.09314 0.1114  0.12286 0.12569 0.12254 0.1127  0.09548 0.06967 0.03804]
  

которая очень похожа на green graph на странице Wiki о бета-биномиальном