Почему моя реализация алгоритма метрополии (mcmc) на python такая медленная?

Question

Я пытаюсь реализовать алгоритм Метрополис (более простая версия алгоритма Метрополис-Гастингс) в Python.

Вот моя реализация:

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

Я ответил на похожий вопрос ранее .Здесь можно использовать многие вещи, которые я упомянул (не вычисление текущей вероятности на каждой итерации, предварительное вычисление случайных инноваций и т. Д.).

Другими улучшениями для вашей реализации было бы не использование списка для храненияваши образцы.Вместо этого вам следует предварительно выделить память для образцов и сохранить их в виде массива.Что-то вроде samples = np.zeros(n_samples) будет более эффективным, чем добавление в список на каждой итерации.

Вы уже упоминали, что пытались повысить эффективность, не записывая записанные пробы.Это хорошая идея.Вы также можете сделать аналогичный трюк с прореживанием, записывая только каждый m-й семпл, так как в любом случае вы отбрасываете их в своем выражении возврата с 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

Если мы посмотрим на производительность, то получим, что эта реализация равна 2xбыстрее , чем оригинал, что неплохо для нескольких незначительных изменений.

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)