Есть много оптимизаций, которые вы можете внести в этот код, прежде чем задуматься о numba et.и др.(Мне удалось ускорить этот код в 25 раз только благодаря умной реализации алгоритма)
Во-первых, в вашей реализации алгоритма Метрополис - Гастингс есть ошибка.Вам нужно сохранять каждую итерацию схемы, независимо от того, движется ваша цепочка или нет.То есть вам нужно удалить posterior = posterior[np.where(posterior > 0)]
из вашего кода и в конце каждого цикла иметь posterior[t] = x_t
.
Во-вторых, этот пример кажется странным.Как правило, с такими типами проблем вывода мы пытаемся вывести параметры распределения, учитывая некоторые наблюдения.Здесь, однако, параметры распределения известны, и вместо этого вы выбираете наблюдения?Как бы то ни было, независимо от этого, я с удовольствием раскрою ваш пример и покажу, как его ускорить.
Ускорение
Для начала удалите все, что не зависитна значение t
из основного цикла for
.Начните с удаления генерации случайного блуждания из цикла 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]
Конечно, также можно переместить случайное поколение u
из цикла 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:
Другая проблема заключается в том, что вы вычисляете текущую апостериорную p2
в каждом цикле, что не является необходимым.В каждом цикле вам нужно рассчитать предложенный апостериорный p1
, и когда предложение будет принято, вы можете обновить 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
Еще одно очень незначительное улучшение заключается в том, что вы заметили, что оператор elif
в вашем коде ничего не делает и поэтому может быть удален.
Сложив все это и используя более разумные имена переменных, я придумал:
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)
Это увеличение скорости на 25x
ESS
Стоит отметить, что грубая скорость - это еще не все, когда речь идет об алгоритмах MCMC.На самом деле, вас интересует количество независимых ( ish ) розыгрышей, которые вы можете сделать сзади в секунду.Как правило, это оценивается с помощью ESS (эффективный размер выборки) .Вы можете повысить эффективность своего алгоритма (и, следовательно, увеличить количество эффективных выборок, извлекаемых в секунду), настроив случайное блуждание.
Для этого обычно выполняется первоначальный пробный запуск, то есть samples = my_get_samples(1000)
.Из этого вывода вычислите sigma = 2.38**2 * np.var(samples)
.Затем это значение следует использовать для настройки случайного блуждания в вашей схеме как innov = norm.rvs(size=n, scale=sigma)
.Казалось бы, произвольное вхождение 2.38 ^ 2 происходит из:
Слабая сходимость и оптимальное масштабирование алгоритмов случайного блуждания Метрополиса (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()
Вы можете сразу увидеть улучшения, внесенные настройкой для повышения эффективности нашего алгоритма.Помните, что оба прогона были сделаны для одинакового количества итераций.![enter image description here](https://i.stack.imgur.com/gf2hE.png)
Заключительные мысли
Если ваш алгоритм требует очень много времени для схождения или если ваши выборки имеют большое количество автокорреляции, я бы подумал посмотретьв Cython
, чтобы выжать дальнейшую оптимизацию скорости.
Я также рекомендовал бы проверить проект PyStan
.Требуется немного привыкнуть, но его алгоритм NUTS HMC, вероятно, превзойдет любой алгоритм Метрополис - Гастингс, который вы можете написать вручную.