对于随机模拟,我需要绘制很多随机数,这是β二项分布的。
目前我以这种方式实现它(使用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毫秒,这在我的情况下主导整个模拟速度。限制元素是β函数(或其中的伽马函数)的评估。
有谁知道如何加快采样速度?
好吧,这里是使用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页面中的绿色图非常相似