【发布时间】：2019-02-19 10:12:30
【问题描述】：
我有一些使用 MCMC 对后验分布进行采样的代码，特别是 Metropolis Hastings。我使用 scipy 生成随机样本：
import numpy as np
from scipy import stats
def get_samples(n):
"""
Generate and return a randomly sampled posterior.
For simplicity, Prior is fixed as Beta(a=2,b=5), Likelihood is fixed as Normal(0,2)
:type n: int
:param n: number of iterations
:rtype: numpy.ndarray
"""
x_t = stats.uniform(0,1).rvs() # initial value
posterior = np.zeros((n,))
for t in range(n):
x_prime = stats.norm(loc=x_t).rvs() # candidate
p1 = stats.beta(a=2,b=5).pdf(x_prime)*stats.norm(loc=0,scale=2).pdf(x_prime) # prior * likelihood
p2 = stats.beta(a=2,b=5).pdf(x_t)*stats.norm(loc=0,scale=2).pdf(x_t) # prior * likelihood
alpha = p1/p2 # ratio
u = stats.uniform(0,1).rvs() # random uniform
if u <= alpha:
x_t = x_prime # accept
posterior[t] = x_t
elif u > alpha:
x_t = x_t # reject
posterior = posterior[np.where(posterior > 0)] # get rid of initial zeros that don't contribute to distribution
return posterior
一般来说，我会尽量避免在 python 中使用显式的 for 循环——我会尝试使用纯 numpy 生成所有内容。然而，对于该算法的情况，带有 if 语句的 for 循环是不可避免的。因此，代码很慢。当我分析我的代码时，大部分时间都在 for 循环中（显然），更具体地说，最慢的部分是生成随机数； stats.beta().pdf() 和 stats.norm().pdf()。
有时我使用numba 来加速我的矩阵运算代码。虽然 numba 与一些 numpy 操作兼容，但生成随机数并不是其中之一。 Numba 有一个cuda rng，但这仅限于正态分布和均匀分布。
我的问题是，有没有办法使用某种与 numba 兼容的各种分布的随机抽样来显着加快上述代码的速度？
我们不必将自己局限于 numba，但它是我所知道的唯一易于使用的优化器。更一般地说，我正在寻找在python中在for循环中加速各种分布（beta、gamma、poisson）随机抽样的方法。
【问题讨论】：
-
Numba 支持很多随机分布。 numba.pydata.org/numba-doc/dev/reference/…我想你必须自己实现的主要是pdf方法（github.com/scipy/scipy/blob/v1.2.1/scipy/stats/…）
-
stats.beta().pdf()和stats.norm().pdf()不会生成随机数。他们正在评估概率分布的密度。 -
顺便说一句，如果您查看 scipy 源代码，您会发现对 scipy.special 函数（xlog1py、betaln）的调用。这是使它们在 Numba 中可用的简单方法github.com/numba/numba/issues/3086
标签： python numpy random numba mcmc