python 中的 beta 二项式分布的有效采样

Question

为了进行随机模拟，我需要绘制大量呈 beta 二项分布的随机数。

目前我是这样实现的（使用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

问题是，抽取一个随机数大约需要 10 分钟。 0.5ms，在我的例子中，这决定了整个模拟速度。限制因素是对 beta 函数（或其中的 gamma 函数）的评估。

有没有人知道如何加快采样速度？

Answer 1

好吧，这是工作代码，经过轻微测试，似乎更快，使用 Beta-Binomial 的复合分布属性。

我们从 beta 中采样 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]

这与 Beta-Binomial 的 Wiki 页面中的绿色图表非常相似

python 中的 beta 二项式分布的有效采样

efficient sampling from beta-binomial distribution in python

python

random

distribution

stochastic

beta-distribution