加速Metropolis - Python中的黑斯廷斯

Question

我有一些代码使用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)的随机抽样.

Answer 1

在开始考虑numba等之前,您可以对此代码进行大量优化.人.(我只是通过智能算法的实现,设法在这段代码上加速了25倍)

首先,您实施Metropolis - Hastings算法时出错.无论您的链是否移动,您都需要保留方案的每次迭代.也就是说,您需要posterior = posterior[np.where(posterior > 0)]从代码中删除并在每个循环结束时都有posterior[t] = x_t.

其次,这个例子看起来很奇怪.通常,对于这些类型的推理问题,我们希望在给出一些观察结果的情况下推断出分布的参数.但是,在这里,分布的参数是已知的,而你是在采样观察？无论如何,无论如何,无论如何,我很高兴能够推出你的例子并向你展示如何加快速度.

加速

首先,您可以从for循环中删除随机游走创新的生成:

    x_t = stats.uniform(0,1).rvs()
    innov = stats.norm(loc=0).rvs(size=n)
    for t in range(n):
        x_prime = x_t + innov[t]

当然也可以t从for循环中移动随机生成:

    x_t = stats.uniform(0,1).rvs()
    innov = stats.norm(loc=0).rvs(size=n)

    u = np.random.uniform(size=n)
    for t in range(n):
        x_prime = x_t + innov[t]
        ...
        if u[t] <= alpha:

另一个问题是你for在每个循环中计算当前后验,这是不必要的.在每个循环中,您需要计算建议的后验u,并且当提议被接受时,您可以更新p2为等于p1:

    x_t = stats.uniform(0,1).rvs()
    innov = stats.norm(loc=0).rvs(size=n)

    u = np.random.uniform(size=n)

    p2 = stats.beta(a=2,b=5).pdf(x_t)*stats.norm(loc=0,scale=2).pdf(x_t)
    for t in range(n):
        x_prime = x_t + innov[t]

        p1 = stats.beta(a=2,b=5).pdf(x_prime)*stats.norm(loc=0,scale=2).pdf(x_prime)
        ...
        if u[t] <= alpha:
            x_t = x_prime # accept
            p2 = p1

        posterior[t] = x_t

一个非常小的改进可能是将scipy stats函数直接导入名称空间:

from scipy.stats import norm, beta

另一个非常小的改进是注意到p2代码中的语句没有做任何事情,因此可以删除.

完全放下这个并使用更明智的变量名称我想出了:

from scipy.stats import norm, beta
import numpy as np

def my_get_samples(n, sigma=1):

    x_cur = np.random.uniform()
    innov = norm.rvs(size=n, scale=sigma)
    u = np.random.uniform(size=n)

    post_cur = beta.pdf(x_cur, a=2, b=5) * norm.pdf(x_cur, loc=0, scale=2)

    posterior = np.zeros(n)
    for t in range(n):
        x_prop = x_cur + innov[t]

        post_prop = beta.pdf(x_prop, a=2, b=5) * norm.pdf(x_prop, loc=0, scale=2)
        alpha = post_prop / post_cur
        if u[t] <= alpha:
            x_cur = x_prop
            post_cur = post_prop

        posterior[t] = x_cur

    return posterior

现在,为了速度比较:

%timeit get_samples(1000)
3.19 s ± 5.28 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit my_get_samples(1000)
127 ms ± 484 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

这是一个25倍的加速

ESS

值得注意的是,当谈到MCMC算法时,粗暴的速度并不是一切.真的,你感兴趣的是你可以从后面每秒进行的独立(ish)绘制的数量.通常,这是使用ESS(有效样本量)评估的.通过调整随机游走,您可以提高算法的效率(从而增加每秒绘制的有效样本).

为此,通常进行初始试运行,即p1.从此输出计算elif.然后应该使用此值来调整方案中的随机游走samples = my_get_samples(1000).似乎任意出现的2.38 ^ 2有它的起源于:

随机游走Metropolis算法的弱收敛和最优缩放(1997).A. Gelman,WR Gilks​​和GO Roberts.

为了说明改进,调整可以让我们运行两次算法,一次调整,另一次调整,都使用10000次迭代:

x = my_get_samples(10000)
y = my_get_samples(10000, sigma=0.12)

fig, ax = plt.subplots(1, 2)
ax[0].hist(x, density=True, bins=25, label='Untuned algorithm', color='C0')
ax[1].hist(y, density=True, bins=25, label='Tuned algorithm', color='C1')
ax[0].set_ylabel('density')
ax[0].set_xlabel('x'), ax[1].set_xlabel('x')
fig.legend()

您可以立即看到调整对我们的算法效率的改进.请记住,两次运行都是针对相同的迭代次数进行的.

最后的想法

如果你的算法花费很长时间来收敛,或者你的样本有大量的自相关,我会考虑sigma = 2.38**2 * np.var(samples)进一步挤出速度优化.

我还建议检查innov = norm.rvs(size=n, scale=sigma)项目.它需要一些习惯,但它的NUTS HMC算法可能会胜过你可以手工编写的任何Metropolis - Hastings算法.