Монте-Карло для генерации точек из распределений (сэмплирование)

Формулировка задачи

Пусть функция \(P(\boldsymbol{r}) \geq 0\) и \(\int d^D \boldsymbol{r} P(\boldsymbol{r})=1\). Тогда её можно трактовать как некоторую плотность распределения.

Хотим генерировать точки \(\boldsymbol{r}_i\) из этого распределения \(P(\boldsymbol{r})\).

При этом считаем, что у нас уже есть способ генерации простых распределений (равномерное, нормальное).

Алгоритм Метрополиса-Хастингса (Markov chain Monte Carlo)

Данный алгоритм основан на случайном блуждании и по сути своей является марковской цепью.

Совсем простыми словами, будем блуждать по области определения, при этом на каждом шаге с бОльшей вероятностью идти в сторону увеличения плотности вероятности. Рекомендую перед прочтением следующего посмотреть вывод кода после параграфа.

Введём \(q\left(x^{\prime} , \; x\right)\) - эмпирическая плотность вероятности перехода из точки \(x\) в точку \(x^{\prime}\). Она не обязано как-то соответствовать целевому распределению \(P(\boldsymbol{r})\); главное условие на \(q\left(x^{\prime} , \; x\right)\) - локальность и умение генерировать из неё точки. Можно думать, что оно вносит консерватизм в систему и не даёт шагам быть слишком большими.

Естественно, если каким-то образом удалось теоретически получить \(q\left(x^{\prime} , \; x\right)\) как истинную вероятность перехода в процессе случайного движения, алгоритм будет работать намного быстрее, но это вовсе необязательно - как правило, хорошо работает многомерный гаусс с маленькой дисперсией.

Сам алгоритм Метрополиса-Хастингса выглядит следующим образом.

1. Выбрать \(x^{\prime}\) по распределению \(q\left(x^{\prime} , \; x^{(t)}\right)\).

Это как бы кандидат на возможное блуждание. К примеру, можно выбирать \(x_{t+1}=x_t+\varepsilon_t\) где \(\varepsilon_t \sim N(0,1)\).

2. Вычислить

\[ a=\frac{P\left(x^{\prime}\right)}{P\left(x^{(t)}\right)} \frac{q\left(x^{(t)} , \; x^{\prime}\right)}{q\left(x^{\prime}, \; x^{(t)}\right)} \]

(не надо пугаться, \(q\left(x^{(t)} , x^{\prime}\right) / q\left(x^{\prime} , x^{(t)}\right)\) для симметричных распределений, например гауссового, равно 1, это просто поправка на асимметрию);

3. С вероятностью а \(\min \left(1, \; a\right)\) генерируем новую точку из распределения \(x^{(t+1)}:=x^{\prime}\), иначе остаёмся на этой же точке \(x^{(t+1)}:=x^{(t)}\).

Итого, результирующий набор \(\{x^{(t)}\}_{t \in T}\) - даст нам произвольную выборку, отвечающую распределению \(P(x)\).

Важные детали:

• Из-за «вольности» в выборке \(q\left(x^{\prime} , \; x^{(t)}\right)\) может возникать корреляция нескольких последующих шагов, что вылезает за применимость цепей Маркова. На практике берут каждый \(n\)-ый шаг в \(\{x^{(t)}\}_{t \in T}\) так, чтобы результирующая выборка отвечала нужным свойствам.

• Начало алгоритма в точке \(x^{(0)}\) нужно стараться выбрать близко примерно в середине \(P(x)\), т.к. неудачный выбор начальной точки может быть сильноскореллированным. На практике первые, скажем, 1000 шагов всегда отбрасывают (burn-in) чтобы точно не словить этот эффект.

• Есть эмпирическое правило - доля отвергнутых кандидатов на блуждание должна быть примерно 0.5 для хорошей скорости. Далеко не всегда работает.

Реализация на python

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.path as path
import scipy.stats as stats

delta = 0.025
X, Y = np.meshgrid(np.arange(-4.5, 2.0, delta), np.arange(-3.5, 3.5, delta))

z1 = stats.multivariate_normal([0, 0], [[1.0, 0], [0, 1.0]])
z2 = stats.multivariate_normal([-2, -2], [[1.5, 0], [0, 0.5]])

def z(x):
    return 0.4 * z1.pdf(x) + 0.6 * z2.pdf(x)

Q = stats.multivariate_normal([0, 0], [[0.05, 0], [0, 0.05]])
r = [0, 0]
samples = [r]
for i in range(0, 1000):
    rq = Q.rvs() + r
    a = z(rq) / z(r)
    if np.random.binomial(1, min(a, 1), 1)[0] == 1:
        r = rq
        samples.append(r)

codes = np.ones(len(samples), int) * path.Path.LINETO
codes[0] = path.Path.MOVETO

p = path.Path(samples, codes)

fig, ax = plt.subplots()

# Add contour lines
Z1 = 0.4 * z1.pdf(np.dstack((X, Y)))
Z2 = 0.6 * z2.pdf(np.dstack((X, Y)))
Z = Z1 + Z2
contour = ax.contour(X, Y, Z, colors='black', linewidths=0.5)

# Add the path
patch = patches.PathPatch(p, facecolor='none', lw=0.5, edgecolor='gray')
ax.add_patch(patch)

ax.set_xlim(-4.5, 2.0)
ax.set_ylim(-3.5, 3.5)
plt.show()

Семплирование по Гиббсу

Пусть размерность пространства \(X\) очень-очень большая. Генерация произвольного блуждания в такой системе может быть очень трудоёмкой (если вообще не невозможной на компьютере, в случае бесконечномерных пространств). Как выпутаться?

Пусть имеем совместное распределение \(P\left(x_1, \ldots, x_D\right)\) для \(D\) случайных величин, причём \(D\) может быть очень большим. Пусть на шаге \(t\) мы уже выбрали какое-то значение \(x_i^t\). На каждом шаге делаются следующие действия:

1. Выбирается индекс \(i:(1 \leq i \leq D)\). Можно идти последовательно по индексам - так сходимость гарантируется.

2. Выбираем \(x_i^{t+1}\) по одномерному распределению \(p\left(x_i \mid x_1^t, \ldots, x_{i-1}^t, x_{i+1}^t, \ldots, x_d^t\right)\) алгоритмом Метрополиса-Хастингса, а для остальных индексов значение не меняется: \(x_j^{t+1}=x_j^t(\mathrm{j} \neq \mathrm{i})\).

Сэмплирование по Гиббсу - это частный случай алгоритма Метрополиса для распределений \(q\left(x^{\prime} , x\right)=p\left(x_i^{\prime} \mid{x}_{-i}\right)\), и вероятность принятия каждого сэмпла получается всегда равна 1 (можете это доказать в качестве упражнения). Поэтому сэмплирование по Гиббсу сходится, и, так как это такое же случайное блуждание по сути, верна та же квадратичная оценка.

Естественно, все ремарки алгоритма Метрополиса-Хастингса насчёт корелированности остаются в силе.

От общего случая блуждания по Метрополису-Гастингсу это отличается тем, что теперь мы двигаемся на каждом шаге вдоль одной из координатных осей; вот соответствующая картинка:

Пример. Генерация гауссова распределения из равномерного

import numpy as np
import matplotlib.pyplot as plt

def proba_density(x):
    
    return 1/np.sqrt(np.pi)*np.exp(-x**2)

x = 0 # Поставим очень плохое начальное приближение
B = 0.1

N = 1000000

xs = np.zeros(N)
xs[0] = x

rejected = 0
for i in range(1, N):
    x_cand = xs[i-1] + np.random.uniform(-B, B)
    
    if proba_density(x_cand)/proba_density(xs[i-1]) >= 1:
        xs[i] = x_cand
        continue
    else:
        if np.random.uniform(0, 1) < proba_density(x_cand)/proba_density(xs[i-1]): # q = 1 т.к. симметричное распределение
            xs[i] = x_cand
            continue
        else:
            xs[i] = xs[i-1]
            rejected  += 1
            continue
plt.hist(xs, bins = 40, density=True)
print(f"Доля отвергнутых кандидатов: {rejected/N}")

Доля отвергнутых кандидатов: 0.02819

eps = 0.5*1e-2
i_cut = 900000

plt.plot(np.linspace(-3, 3, 200), proba_density(np.linspace(-3, 3, 200)))
_ = plt.hist(xs[i_cut::3], bins = 40, density = True)
xs

array([ 0.        ,  0.03577155, -0.04432745, ..., -1.05114796,
       -1.03528899, -1.09977569])