эффективная выборка из бета-биномиального распределения в Python - PullRequest
Новая
       75

эффективная выборка из бета-биномиального распределения в Python

2 голосов
/ 12 марта 2019

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

На данный момент я реализовал это так (используя 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 мс, что в моем случае доминирует над всей скоростью симуляции. Ограничивающим элементом является оценка бета-функций (или гамма-функций внутри них).

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

1 Ответ

0 голосов
/ 16 марта 2019

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

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

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]

, что очень похоже на зеленый график на странице Wiki в бета-биномиальном виде

...