Есть много оптимизаций, которые вы можете внести в этот код, прежде чем задуматься о 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()
Вы можете сразу увидеть улучшения, внесенные настройкой для повышения эффективности нашего алгоритма.Помните, что оба прогона были сделаны для одинакового количества итераций.
Заключительные мысли
Если ваш алгоритм требует очень много времени для схождения или если ваши выборки имеют большое количество автокорреляции, я бы подумал посмотретьв Cython, чтобы выжать дальнейшую оптимизацию скорости.
Я также рекомендовал бы проверить проект PyStan.Требуется немного привыкнуть, но его алгоритм NUTS HMC, вероятно, превзойдет любой алгоритм Метрополис - Гастингс, который вы можете написать вручную.