python - 为什么我的大都会算法(mcmc)的Python实现这么慢？

Question

我正在尝试实现 Metropolis Python 中的算法(Metropolis-Hastings 算法的简单版本)。

这是我的实现:

def Metropolis_Gaussian(p, z0, sigma, n_samples=100, burn_in=0, m=1):
    """
    Metropolis Algorithm using a Gaussian proposal distribution.
    p: distribution that we want to sample from (can be unnormalized)
    z0: Initial sample
    sigma: standard deviation of the proposal normal distribution.
    n_samples: number of final samples that we want to obtain.
    burn_in: number of initial samples to discard.
    m: this number is used to take every mth sample at the end
    """
    # List of samples, check feasibility of first sample and set z to first sample
    sample_list = [z0]
    _ = p(z0) 
    z = z0
    # set a counter of samples for burn-in
    n_sampled = 0

    while len(sample_list[::m]) < n_samples:
        # Sample a candidate from Normal(mu, sigma),  draw a uniform sample, find acceptance probability
        cand = np.random.normal(loc=z, scale=sigma)
        u = np.random.rand()
        try:
            prob = min(1, p(cand) / p(z))
        except (OverflowError, ValueError) as error:
            continue
        n_sampled += 1

        if prob > u:
            z = cand  # accept and make candidate the new sample

        # do not add burn-in samples
        if n_sampled > burn_in:
            sample_list.append(z)

    # Finally want to take every Mth sample in order to achieve independence
    return np.array(sample_list)[::m]

当我尝试将我的算法应用于指数函数时，只需要很少的时间。然而，当我在t-分布上尝试它时，考虑到它没有进行那么多计算，它需要很长时间。这是复制我的代码的方法:

t_samples = Metropolis_Gaussian(pdf_t, 3, 1, 1000, 1000, m=100)
plt.hist(t_samples, density=True, bins=15, label='histogram of samples')
x = np.linspace(min(t_samples), max(t_samples), 100)
plt.plot(x, pdf_t(x), label='t pdf')
plt.xlim(min(t_samples), max(t_samples))
plt.title("Sampling t distribution via Metropolis")
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.legend()

这段代码需要相当长的时间才能运行，我不知道为什么。在我的 Metropolis_Gaussian 代码中，我试图通过以下方式提高效率

1. 不将重复样本添加到列表中
2. 不记录烧机样本

<小时/>

函数pdf_t定义如下

from scipy.stats import t
def pdf_t(x, df=10):
    return t.pdf(x, df=df)

Answer 1

我回答了similar question previously 。我在那里提到的许多东西(不是在每次迭代时计算当前的可能性，预先计算随机创新等)都可以在这里使用。

实现的其他改进是不使用列表来存储样本。相反，您应该为样本预先分配内存并将它们存储为数组。像 samples = np.zeros(n_samples) 这样的东西比在每次迭代时附加到列表更有效。

您已经提到您试图通过不记录老化样本来提高效率。这是一个好主意。您还可以通过仅记录每个第 m 个样本来对细化执行类似的技巧，因为您无论如何都会在 return 语句中使用 np.array(sample_list)[::m] 丢弃这些样本。您可以通过更改来做到这一点:

   # do not add burn-in samples
    if n_sampled > burn_in:
        sample_list.append(z)

至

    # Only keep iterations after burn-in and for every m-th iteration
    if n_sampled > burn_in and n_sampled % m == 0:
        samples[(n_sampled - burn_in) // m] = z

还值得注意的是，您不需要计算 min(1, p(cand)/p(z))，只需计算 p(cand)/p(z) 。我意识到，形式上 min 是必要的(以确保概率限制在 0 和 1 之间)。但是，在计算上，我们不需要最小值，因为如果 p(cand)/p(z) > 1 那么 p(cand)/p(z) 是总是大于u。

将这些放在一起，并预先计算随机创新、接受概率u，并且仅在确实需要时才计算可能性，我想出了:

def my_Metropolis_Gaussian(p, z0, sigma, n_samples=100, burn_in=0, m=1):
    # Pre-allocate memory for samples (much more efficient than using append)
    samples = np.zeros(n_samples)

    # Store initial value
    samples[0] = z0
    z = z0
    # Compute the current likelihood
    l_cur = p(z)

    # Counter
    iter = 0
    # Total number of iterations to make to achieve desired number of samples
    iters = (n_samples * m) + burn_in

    # Sample outside the for loop
    innov = np.random.normal(loc=0, scale=sigma, size=iters)
    u = np.random.rand(iters)

    while iter < iters:
        # Random walk innovation on z
        cand = z + innov[iter]

        # Compute candidate likelihood
        l_cand = p(cand)

        # Accept or reject candidate
        if l_cand / l_cur > u[iter]:
            z = cand
            l_cur = l_cand

        # Only keep iterations after burn-in and for every m-th iteration
        if iter > burn_in and iter % m == 0:
            samples[(iter - burn_in) // m] = z

        iter += 1

    return samples

如果我们看一下性能，我们会发现此实现比原始实现快 2 倍，这对于一些小的更改来说还不错。

In [1]: %timeit Metropolis_Gaussian(pdf_t, 3, 1, n_samples=100, burn_in=100, m=10)
205 ms ± 2.16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [2]: %timeit my_Metropolis_Gaussian(pdf_t, 3, 1, n_samples=100, burn_in=100, m=10)
102 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)