notes/machine_learning/bayesian_methods.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы"
]
},
"source": [
"## Байесовский подход к дискриминативным моделям\n",
"\n",
"Под априорной (необученной) дискриминативной моделью распределения матрицы «объекты-признаки» $X$, вектора целевой переменной $y$ и параметров модели $\\theta$ будем понимать условную плотность $p(y, \\theta \\, \\vert \\, X) = p(y \\, \\vert \\, X, \\theta) p(\\theta \\, \\vert X) = p(y \\, \\vert \\, X, \\theta) p(\\theta)$, где последнее равенство верно, потому что предполагается, что априорное распределение параметров не зависит от данных. Доступность такой априорной дискриминативной модели будем понимать как то, что оба множителя из правой части (то есть $p(y \\, \\vert \\, X, \\theta)$ и $p(\\theta)$) известны в аналитическом виде.\n",
"\n",
"В байесовском смысле задачей обучения дискриминативной модели является задача нахождения апостериорной условной плотности $p(\\theta \\, \\vert \\, X_\\mathrm{train}, y_\\mathrm{train})$, где $X_\\mathrm{train}$ и $y_\\mathrm{train}$ образуют обучающую выборку. Искомая условная плотность находится по теореме Байеса, где «меняются местами» $\\theta$ и $y_\\mathrm{train}$ и где все вероятности являются условными относительно $X_\\mathrm{train}$ (просто на заметку: в дискриминативном моделировании в отличие от генеративного моделирования не предполагается, что известна $p(X_\\mathrm{train})$):\n",
"$$p(\\theta \\, \\vert \\, X_\\mathrm{train}, y_\\mathrm{train}) = \\frac{p(y_\\mathrm{train} \\, \\vert \\, X_\\mathrm{train}, \\theta) p(\\theta)}{\\int p(y_\\mathrm{train} \\, \\vert \\, X_\\mathrm{train}, \\theta) p(\\theta) d\\theta}.$$\n",
"В предположении, что интеграл из знаменателя нижеследующего равенства берётся, задача обучения решена, а если он не берётся, то используются приближённые байесовские методы.\n",
"\n",
"Применение же обученной модели в байесовском смысле является задачей нахождения распределения $p(y_\\mathrm{new} \\, \\vert \\, X_\\mathrm{new}, X_\\mathrm{train}, y_\\mathrm{train})$ вектора целевой переменной $y_\\mathrm{new}$ для известной матрицы «объекты-признаки» $X_\\mathrm{new}$. На стадии применения возникает взвешенный ансамбль моделей с различными параметрами:\n",
"$$p(y_\\mathrm{new} \\, \\vert \\, X_\\mathrm{new}, X_\\mathrm{train}, y_\\mathrm{train}) = \\int p(y_\\mathrm{new} \\, \\vert \\, X_\\mathrm{new}, \\theta) p(\\theta \\, \\vert \\, X_\\mathrm{train}, y_\\mathrm{train}) d\\theta,$$\n",
"где плотность $p(\\theta \\, \\vert \\, X_\\mathrm{train}, y_\\mathrm{train})$ была найдена при обучении, а если бы вместо неё стояла априорная плотность $p(\\theta)$, то слева получилась бы плотность $p(y_\\mathrm{new} \\, \\vert \\, X_\\mathrm{new})$."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы",
"выбор_модели"
]
},
"source": [
"## Принцип наибольшей обоснованности\n",
"\n",
"Зададимся вопросом: если есть множество из $m$ [дискриминативных моделей](__home_url__/notes/Байесовский подход к дискриминативным моделям) $\\{p_i(y \\, \\vert \\, X, \\theta_i) p_i(\\theta_i)\\}_{i=1}^m$, то какую из них выбрать? Один из ответов даёт принцип наибольшей обоснованности (evidence), предлагающий руководствоваться следующим правилом:\n",
"$$i_\\mathrm{opt} = \\arg \\max_i \\int p_i(y_\\mathrm{train} \\, \\vert \\, X_\\mathrm{train}, \\theta) p_i(\\theta) d\\theta,$$\n",
"где за $\\theta$ обозначен вектор параметров из множества, являющегося объединением множеств параметров всех $m$ моделей (если у $i$-й модели нет какого-то параметра из $\\theta$, то можно, например, считать, что в рамках $i$-й модели его распределение задаётся дельта-функцией Дирака).\n",
"\n",
"С философской точки зрения, принцип наибольшей обоснованности является математической формализацией бритвы Оккама: из моделей, способных объяснить имеющиеся данные, выбирается наиболее простая, то есть наименее способная объяснять при другом выборе параметров данные, не похожие на имеющиеся. Если же говорить более предметно, то можно сказать, что принцип наибольшей обоснованности ищет баланс между двумя желательными, но не всегда совместимыми друг с другом свойствами модели. Эти свойства таковы:\n",
"* способность хорошо воспроизводить зависимость, присутствующую в обучающей выборке: чем выше в среднем $p_i(y_\\mathrm{train} \\, \\vert \\, X_\\mathrm{train}, \\theta)$, тем выше обоснованность;\n",
"* неспособность к переподгонке (overfitting): в общем случае, в зависимости от $\\theta$ модель может воспроизводить различные зависимости и чем вероятнее $\\theta$, при которых получается какая-то зависимость, не согласующаяся с данными, тем ниже обоснованность. Отметим, что усреднение происходит относительно априорного распределения $p(\\theta)$, так что даже если есть значение вектора $\\theta$, для которого высока $p_i(y_\\mathrm{train} \\, \\vert \\, X_\\mathrm{train}, \\theta)$ и которое поэтому стало бы наиболее вероятным после обучения, вклад такого $\\theta$ будет внесён лишь с весом $p(\\theta)$, а не с весом $p(\\theta \\, \\vert \\, X_\\mathrm{train}, y_\\mathrm{train})$.\n",
"\n",
"Считая, что модели отличаются только гиперпараметрами, при помощи принципа наибольшей обоснованности можно подбирать гиперпараметры без кросс-валидации. Это особенно актуально тогда, когда гиперпараметров много, из-за чего поиск по хоть сколько-нибудь мелкой сетке становится вычислительно неподъёмным.\n",
"\n",
"Более подробно про принцип наибольшей обоснованности можно прочесть в главе 28 книги [MacKay, 2003](http://www.inference.org.uk/mackay/itprnn/ps/343.355.pdf)."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы",
"статистический_вывод",
"EM"
]
},
"source": [
"## Общий взгляд на EM-алгоритм\n",
"\n",
"#### Введение\n",
"\n",
"В задаче кластеризации EM-алгоритм может использоваться для разделения смеси гауссовских распределений. Он имеет ряд преимуществ перед методом K-средних. Во-первых, вместо бинарных индикаторов принадлежности к кластеру вводятся нечёткие степени принадлежности к кластеру. Во-вторых, ковариационные матрицы гауссиан не обязаны быть единичными (как это неявно предполагается при измерении евклидова расстояния до центров кластеров в алгоритме K-средних)— для отсутствия вычислительной неустойчивости достаточно, чтобы они были диагональными, хотя есть и версии с более слабыми ограничениями.\n",
"\n",
"На самом же деле, EM-алгоритм имеет более общий вид, чем алгоритм для разделения смеси гауссовских распределений. Это алгоритм, позволяющий делать две вещи:\n",
"* оценивать параметры распределений тогда, когда часть переменных относится к ненаблюдаемым;\n",
"* решать задачу статистического вывода (statistical inference), то есть задачу оценивания значений ненаблюдаемых признаков объекта по его значениям наблюдаемых признаков.\n",
"\n",
"Более того, EM-алгоритм может применяться и к выборкам, где векторы объектов не являются независимыми одинаково распределёнными: например, он может применяться к последовательностям (таким как марковские цепи).\n",
"\n",
"#### Постановка задачи\n",
"\n",
"Формализуем постановку задачи так.\n",
"\n",
"Пусть есть следующие сущности:\n",
"* $X$ — $(l \\times n)$-матрица наблюдаемых переменных, где $l$ — количество объектов в выборке, а $n$ — количество наблюдаемых признаков;\n",
"* $T$ — $(l \\times h)$-матрица скрытых переменных, где $h$ — количество скрытых признаков, а $l$ то же самое, что и в матрице $X$. Для получения полных признаковых описаний объектов выборки матрицу $X$ нужно соединить с матрицей $T$ боковой конкатенацией;\n",
"* $\\theta$ — неизвестный вектор параметров вероятностного распределения, породившего выборку; этот вектор требуется оценить по данным; в обычной постановке это детерминированный вектор, в байесовской же постановке (рассматриваемой лишь в самом конце заметки) это случайный вектор со своим собственным распределением.\n",
"\n",
"Эти три сущности связаны между собой через функцию совместной условной плотности $p(X, T \\, \\vert \\, \\theta)$. В рассматриваемой задаче эта функция дана в аналитическом виде для любой длины $l$ матриц $X$ и $T$.\n",
"\n",
"Тут уместно пояснить, почему выше используются матрицы $X$ и $T$ вместо векторов-строк $x$ и $t$, соответствующих какому-либо одному объекту. Дело в том, что различные строки матриц $X$ или $T$ могут не быть независимыми друг от друга, так что в общем случае функцию плотности нужно брать не от отдельных объектов выборки, а от всей выборки сразу. И, на самом деле, постановку можно ещё более обобщить, определив $X$ как кортеж (tuple) наблюдаемых переменных, $T$ — как кортеж скрытых переменных, а $\\theta$ — как кортеж параметров модели.\n",
"\n",
"Задача заключается в том, чтобы найти оценку для $\\theta$ методом максимизации неполного правдоподобия, то есть правдоподобия только матрицы наблюдаемых переменных $X$, а не совместного правдоподобия $X$ и $T$: \n",
"$$\\hat{\\theta} = \\arg\\max_\\theta p(X \\, \\vert \\, \\theta) = \\arg\\max_\\theta \\int p(X, T \\, \\vert \\, \\theta) dT.$$\n",
"Неполное правдоподобие используется, потому что значение матрицы $T$ не дано, а дано только значение матрицы $X$.\n",
"\n",
"Имея оценку $\\hat{\\theta}$, можно по матрице с наблюдаемой частью признакового описания объектов $X$ выводить матрицу ненаблюдаемой части признакового описания этих объектов $T$ следующим образом:\n",
"$$\\hat{T} = \\arg\\max_T p(X, T \\, \\vert \\, \\hat{\\theta}).$$\n",
"Это точечная оценка, а из описания E-шага станет понятно, как получать оценки распределения матрицы $T$.\n",
"\n",
"#### Пример частной постановки\n",
"\n",
"Задача разделения смеси гауссовских распределений является частным случаем описанной задачи, получающимся из неё, если положить, что:\n",
"* $X$ — $(l \\times n)$-матрица наблюдаемых признаков объектов;\n",
"* $T$ — $(l \\times 1)$-матрица (т.е. вектор-столбец) с дискретной скрытой переменной, задающей идентификатор (номер) компоненты смеси, откуда был порождён объект;\n",
"* $\\theta$ — неизвестный вектор, составленный из:\n",
" - вероятностей компонент смеси $p(i)$, где $i$ — скалярный идентификатор компоненты смеси,\n",
" - центров гауссовских распределений $(\\mu_i)_{i=1}^m$, где $m$ — количество компонент смеси,\n",
" - ковариационных матриц гауссовских распределений $(S_i)_{i=1}^m$, развёрнутых из двумерного в одномерный вид, чтобы стать компонентами вектора $\\theta$.\n",
"\n",
"Чтобы показать, что задача разделения смеси гауссовских распределений является частным случаем поставленной задачи, необходимо и достаточно убедиться, что $p(X, T \\, \\vert \\, \\theta)$ дана в аналитическом виде. Поскольку в задаче разделения смеси гауссовских распределений объекты являются независимыми одинаково распределёнными, верно, что:\n",
"$$p(X, T \\, \\vert \\, \\theta) = \\prod_{i=1}^l p(x_{i\\cdot}, t_{i\\cdot} \\, \\vert \\, \\theta),$$\n",
"где $x_{i\\cdot}$ обозначает $i$-ю строку матрицы $X$, а $t_{i\\cdot}$ обозначает $i$-ю строку матрицы $T$. Для векторов-строк $x_{i\\cdot}$ и $t_{i\\cdot}$ (т.е. матриц размера $1 \\times n$ и $1 \\times 1$ соответственно) плотность имеет следующий вид: \n",
"$$p(x_{i\\cdot}, t_{i\\cdot} \\, \\vert \\, \\theta) = p(t_{i\\cdot}) \\mathcal{N}_{\\mu_{t_{i\\cdot}}, S_{t_{i\\cdot}}}(x_{i\\cdot}),$$\n",
"где в $(1 \\times 1)$-матрице $t_{i\\cdot}$ хранится скалярный идентификатор компоненты смеси, $\\mathcal{N}_{\\mu_{t_{i\\cdot}}, S_{t_{i\\cdot}}}(x_{i\\cdot})$ — плотность гауссовского распределения с соответствующими параметрами в точке с координатами $x_{i\\cdot}$, а $p(t_{i\\cdot})$, $\\mu_{t_{i\\cdot}}$ и $S_{t_{i\\cdot}}$ извлечены из $\\theta$. Подстановкой этого равенства в предыдущее и получается аналитический вид для $p(X, T \\, \\vert \\, \\theta)$.\n",
"\n",
"#### Обозначения\n",
"\n",
"Обозначим за $q(T)$ плотность некоторого (не обязательно известного) вероятностного распределения матрицы $T$. В дальнейшем будет считаться, что эта плотность приближает плотность апостериорной вероятности $p(T \\, \\vert \\, X, \\theta)$. Например, в задаче разделения смеси гауссовских распределений $q(T)$ хранит информацию о сформировавшихся к текущей итерации EM-алгоритма нечётких степенях принадлежности к кластерам.\n",
"\n",
"Наконец, для определённости в том, использовать ли знак интеграла или знак суммы, будем считать, что $T$ является матрицей с элементами из $\\mathbb{R}$. Если элементы $T$ дискретны, то интегралы заменятся на суммы, но содержательно ничего не изменится.\n",
"\n",
"#### Нижняя вариационная оценка\n",
"\n",
"Для решения общей задачи предпримем следующие шаги, на каждом из которых будем преобразовывать логарифм правдоподобия наблюдаемой выборки:\n",
"* $\\log p(X \\, \\vert \\, \\theta) = \\int q(T) \\log p(X \\, \\vert \\, \\theta) dT$, потому что левая часть не зависит от $T$, а справа её просто внесли как константу в интеграл, равный 1;\n",
"* правую часть только что выписанного равенства представим в виде $\\int q(T) \\log \\frac{p(X, T \\, \\vert \\, \\theta)}{p(T \\, \\vert \\, X, \\theta)} dT$ — сделать это можно по определению условной вероятности;\n",
"* домножим выражение под логарифмом на единицу, домножив числитель и знаменатель на одно и то же: $\\int q(T) \\log \\frac{p(X, T \\, \\vert \\, \\theta) \\, q(T)}{p(T \\, \\vert \\, X, \\theta) \\, q(T)} dT$;\n",
"* воспользуемся тем, что логарифм произведения равен сумме логарифмов, а множители сгруппируем крест-накрест: $\\int q(T) \\log \\frac{p(X, T \\, \\vert \\, \\theta)}{q(T)} dT + \\int q(T) \\log \\frac{q(T)}{p(T \\, \\vert \\, X, \\theta)} dT$.\n",
"\n",
"Первое слагаемое обозначим за $\\mathcal{L}(q, \\theta)$:\n",
"$$\\mathcal{L}(q, \\theta) = \\int q(T) \\log \\frac{p(X, T \\, \\vert \\, \\theta)}{q(T)} dT.$$\n",
"Отметим, что в правой части присутствует $X$, но он не выписан как аргумент, от которого зависит $\\mathcal{L}(q, \\theta)$. Опустить $X$ можно по следующей причине: это уже заданная константная матрица пронаблюдавшихся значений. Таким образом, $\\mathcal{L}(q, \\theta)$ является функционалом, зависящим от функции плотности $q$ и вектора $\\theta$.\n",
"\n",
"Кроме того заметим, что второе слагаемое в получившемся на последнем шаге выражении равно [дивергенции Кульбака-Лейблера](__home_url__/notes/Дивергенция Кульбака-Лейблера) между некоторым распределением $q(T)$ и распределением $p(T \\, \\vert \\, X, \\theta)$, которое им приближается по предположению, сделанному при описании $q(T)$ (впрочем, порядок аргументов дивергенции Кульбака-Лейблера такой, как если бы, наоборот, распределение $q(T)$ приближалось распределением $p(T \\, \\vert \\, X, \\theta)$):\n",
"$$D_{KL}(q(T) \\, \\Vert \\, p(T \\, \\vert \\, X, \\theta)) = \\int q(T) \\log \\frac{q(T)}{p(T \\, \\vert \\, X, \\theta)} dT.$$\n",
"\n",
"Итак, путём вышеописанных преобразований было получено, что:\n",
"$$\\log p(X \\, \\vert \\, \\theta) = \\mathcal{L}(q, \\theta) + D_{KL}(q(T) \\, \\Vert \\, p(T \\, \\vert \\, X, \\theta)).$$\n",
"Второе слагаемое правой части неотрицательно в силу соответствующего свойства дивергенции Кульбака-Лейблера. Из этих соображений первое слагаемое называют нижней вариационной оценкой для $\\log p(X \\, \\vert \\, \\theta)$.\n",
"\n",
"#### Описание EM-алгоритма\n",
"\n",
"Теперь всё готово, для того чтобы описать EM-алгоритм, пошагово максимизирующий то по $q$, то по $\\theta$ прологарифмированное правдоподобие $\\log p(X \\, \\vert \\, \\theta)$:\n",
"* E-шаг: при фиксированном $\\theta = \\theta^{(n)}$ решается задача $\\mathcal{L}(q, \\theta^{(n)}) \\to \\max_{q \\in Q}$, где $Q$ — некоторое заранее выбранное множество вероятностных распределений. Для получения общего представления о том, как проводится такая максимизация, рассмотрим случай, когда $Q$ — множество всех вероятностных распределений. В этом случае делаются точные E-шаги, а решение оптимизационной задачи получается тогда, когда обнуляется дивергенция Кульбака-Лейблера, поскольку $\\mathcal{L}(q, \\theta^{(n)})$ не может быть больше $\\log p(X \\, \\vert \\, \\theta^{(n)})$, в рамках текущего шага являющегося константой. Значит, решением задачи максимизации является распределение $q^{(n+1)}(T) = p(T \\, \\vert \\, X, \\theta^{(n)}) = \\frac{p(X, T \\, \\vert \\, \\theta^{(n)})}{\\int p(X, T \\, \\vert \\, \\theta^{(n)}) dT}$, где правая часть (с точностью до взятия интеграла из знаменателя) известна, потому что $X$ известен, $\\theta^{(n)}$ известен и по предоположению, сделанному в формулировке задачи, известен аналитический вид $p(X, T \\, \\vert \\, \\theta)$. Так как параметр $\\theta = \\theta^{(n)}$ зафиксирован, можно сказать, что из параметрического семейства распределений $p_{\\theta}(X, T) = p(X, T \\, \\vert \\, \\theta)$ выбрано какое-то одно $p(X, T \\, \\vert \\, \\theta^{(n)})$ и то, что делается на E-шаге, является решением задачи [статистического вывода](__home_url__/notes/Байесовский статистический вывод) для этого распределения;\n",
"* М-шаг: при фиксированном $q = q^{(n)}$ решается задача $\\mathcal{L}(q^{(n)}, \\theta) \\to \\max_{\\theta}$. Поскольку логарифм частного равен разности логарифмов делимого и делителя, а распределение $q^{(n)}$ зафиксировано, задача эквивалентна задаче $\\mathbb{E}_{T \\sim q^{(n)}} \\log p(X, T \\, \\vert \\, \\theta) \\to \\max_{\\theta}$, где матрица $T$ случайна и именно по её распределению берётся математическое ожидание. Вообще говоря, выписанная оптимизационная задача может быть многоэкстремальной, но если $Q$ — экспоненциальный класс распределений, то задача выпуклая (максимизируется вогнутая функция).\n",
"\n",
"#### Различные версии EM-алгоритма\n",
"\n",
"В зависимости от вида ограничения $q \\in Q$ возникают различные версии EM-алгоритма:\n",
"* Точный EM-алгоритм. В нём $Q$ — множество всех вероятностных распределений. Этот пример уже разбирался в описании EM-алгоритма. Осталось лишь заметить, что на E-шаге при вычислении решения требуется взять интеграл $\\int p(X, T \\, \\vert \\, \\theta^{(n)}) dT$, но он может не браться. Проблема может возникнуть, если распределения $p(X \\, \\vert \\, T, \\theta^{(n)})$ и $p(T \\, \\vert \\, \\theta^{(n)})$ не образуют пару сопряжённых друг к другу. Поэтому иногда приходится брать как $Q$ нечто более узкое, чем множество всех вероятностных распределений, и проводить приближённые E-шаги.\n",
"* Жёсткий (crisp) EM-алгоритм. В нём $Q$ — множество всех атомарных распределений, сконцентрированных в одной точке (в данном случае, в одной матрице): $Q = \\{\\delta(T - T_0) \\, \\vert \\, T_0 \\in \\mathbb{R}^{l \\times h}\\}$, где $\\delta$ — дельта-функция Дирака. При таком ограничении E-шаг всегда можно провести: из определения дивергенции Кульбака-Лейблера получится, что решение — это распределение, сконцентрированное в точке моды по $T$ распределения $p(T \\, \\vert \\, X, \\theta^{(n)}))$, то есть в точке $\\arg \\max_{T} p(T \\, \\vert \\, X, \\theta^{(n)}))$. Говоря абстрактно, можно сказать, что в общем случае задача E-шага является задачей функциональной оптимизации (оптимум ищется по вероятностному распределению $q$), а такие задачи чаще всего не могут быть решены средствами современной математики, однако сужением $Q$ до параметрического семейства задачу можно свести к задаче оптимизации по вещественным параметрам (например, все элементы семейства атомарных распределений биективно сопоставляются координатам точки/матрицы, где сконцентрирована вся масса).\n",
"* Вариационный EM-алгоритм. В нём $Q$ — множество распределений, допускающих факторизацию по строкам: $Q = \\{ q(T)\\, \\vert \\, q(T) = \\prod_{j=1}^l q_j(t_{j\\cdot})\\}$, где $q_j(t_{j\\cdot})$ — некоторое распределение уже только одного вектора-строки, соответствующего $j$-му наблюдению (на самом деле, такое предположение приводит к naive mean-field approximation, но можно ослабить ограничение, вместо отдельных строк перейдя к непересекающимся множествам строк, что приведёт к generalized mean-field approximation). Заметим, что дельта-функция от матрицы раскладывается в произведение дельта-функций от своих строк, так что вариационный EM-алгоритм охватывает все варианты, доступные жёсткому EM-алгоритму. На E-шаге решение $q^{(n+1)}$ ищется с использованием вариационного подхода (покоординатного подъёма), и при соблюдении некоторых условий это решение можно найти.\n",
"\n",
"Наконец, отдельно стоит выделить байесовский EM-алгоритм, где $\\theta$ рассматривается как случайный вектор, имеющий своё собственное известное априорное распределение $p(\\theta)$, а оценить требуется апостериорные распределения $p(\\theta \\, \\vert \\, X)$ и $p(T \\, \\vert \\, X, \\theta)$. Введение $p(\\theta)$ влечёт за собой то, что вместо $p(X, T \\, \\vert \\, \\theta)$ в выкладках фигурирует $p(X, T, \\theta) = p(X, T \\, \\vert \\, \\theta) p(\\theta)$. E-шаг при этом не изменится, потому что $p(\\theta^{(n)})$ сократится из выражения для $p(T \\, \\vert \\, X, \\theta^{(n)})$, а вот на M-шаге будет решаться следующая задача:\n",
"$$\\mathbb{E}_{T \\sim q^{(n)}} \\log p(X, T, \\theta) \\to \\max_{\\theta}.$$\n",
"Иными словами, на M-шаге по сравнению с версией без априорного распределения на $\\theta$ к оптимизируемому выражению просто добавится слагаемое $\\log p(\\theta)$.\n",
"\n",
"#### Преимущества EM-алгоритма\n",
"\n",
"Вспомним поставновку задачи, для решения которой был предложен EM-алгоритм и зададимся вопросом, почему нельзя было просто максимизировать неполное правдоподобие так же, как максимизируют обычное правдоподобие. Ответ заключается в следующем:\n",
"* Хотя подынтегральное выражение $p(X, T \\, \\vert \\, \\theta)$ и дано в аналитическом виде, часто бывает так, что сам интеграл $\\int p(X, T \\, \\vert \\, \\theta) dT$ в аналитическом виде не берётся, так что брать градиент не от чего.\n",
"* Даже если неполное правдоподобие удаётся вычислить аналитически, EM-алгоритм, как правило, работает быстрее, чем градиентный спуск. Конечно, если $l$ велико, то используют стохастический градиентный спуск, но существует и стохастический EM-алгоритм, где на каждой итерации используется лишь часть выборки.\n",
"* На выходе также получается информация о $T$, а это может помочь лучше понять данные.\n",
"* Если часть элементов матрицы $X$ пропущена, не нужно применять эвристики, не имеющие под самой математического обоснования (например, такие как замена пропусков на среднее по столбцу). Просто нужно немного изменить постановку, к наблюдаемым переменным отнеся только то, что известно, и включив пропуски в кортеж скрытых переменных."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы",
"статистический_вывод"
]
},
"source": [
"## Байесовский статистический вывод\n",
"\n",
"#### Задача статистического вывода\n",
"\n",
"Пусть имеются следующие сущности:\n",
"* кортеж (tuple) наблюдаемых переменных $X$, индексируемый одним индексом $i$; чаще всего будем считать, что этот кортеж составлен из известной выборки, то есть из матрицы «объекты-признаки», содержащей только наблюдаемые признаки объектов и имеющей размеры $l \\times n$, где $l$ — количество объектов, а $n$ — количество наблюдаемых признаков, в таком случае запись $X_i$ будет обозначать $i$-ю строку этой матрицы;\n",
"* кортеж скрытых переменных $Z$, индексируемый одним индексом $i$; в $Z$ может входить матрица скрытых признаков объектов $T$ размера $l \\times h$, где $h$ — количество скрытых признаков, и тогда запись $Z_i$, где $1 \\le i \\le l$, обозначает $i$-ю строку матрицы $T$, то есть вектор-строку всех скрытых признаков $i$-го объекта; также в $Z$ могут входить неизвестные параметры модели, и тогда запись $Z_i$, где $l < i$, обозначает $(i-l)$-й параметр модели;\n",
"* вероятностная модель, то есть аналитический вид совместного распределения $p(X, Z)$.\n",
"\n",
"Плотность вероятностного распределения дана именно для кортежей $X$ и $Z$, а не для одного произвольного объекта и параметров модели по ряду причин: в том числе по той, что в общем случае объекты могут не быть независимы друг от друга. Например, ими могут быть наблюдения элементов одной и той же последовательности с влиянием прошлого на будущее.\n",
"\n",
"Задача заключается в нахождении оценки $\\hat{Z}$ для неизвестного кортежа $Z$, соответствующего известному кортежу $X$. Также задачу можно расширить до такой: оценить апостериорное распределение $p(Z \\, \\vert \\, X)$, а не только найти точечную оценку $\\hat{Z}$.\n",
"\n",
"#### Пример\n",
"\n",
"Иллюстративным примером заадчи статистического вывода является задача байесовского обучения линейной регрессии, сформулированная в терминах нижеследующих сущностей:\n",
"* $X$ — кортеж из обучающей выборки, то есть из строк матрицы «объекты-признаки», где помимо собственно признаков есть и целевая переменная;\n",
"* $Z$ — кортеж $(w, \\alpha, \\beta)$, где $w$ — $(n \\times 1)$-вектор весов линейной модели, $\\alpha$ — скалярный гиперпараметр регуляризации, а $\\beta$ — единица, делённая на дисперсию шума, присутствующего в наблюдениях;\n",
"* вероятностная модель имеет вид $p(X, Z) = p(X \\, \\vert \\, Z) p(Z)$, где:\n",
" - $p(X \\, \\vert \\, Z) = \\prod_{i=1}^l \\mathcal{N}_{x_i w, \\beta^{-1}}(y_i)$, где $\\mathcal{N}_{x_i w, \\beta^{-1}}(y_i)$ — плотность в точке $y_i$ нормального распределения со средним $x_i w$ и дисперсией $\\beta^{-1}$, $x_i$ — $(1 \\times n)$-вектор только признаков $i$-го объекта, $y_i$ — значение целевой переменной на $i$-м объекте;\n",
" - $p(Z) = p(w \\, \\vert \\, \\alpha) p(\\alpha) p(\\beta)$,\n",
" - $p(w \\, \\vert \\, \\alpha) = \\mathcal{N}_{0, \\alpha^{-1}I}(w)$, где $I$ — единичная матрица размера $n \\times n$,\n",
" - $p(\\alpha) = \\Gamma_{a, b}(\\alpha)$, $p(\\beta) = \\Gamma_{c, d}(\\beta)$, где $\\Gamma$ — плотность гамма-распределения, а $a$, $b$, $c$ и $d$ — константные гиперпараметры априорных распределений.\n",
"\n",
"#### Байесовский подход к задаче статистического вывода\n",
"\n",
"При байесовском подходе предполагается, что $p(X, Z) = p(X \\, \\vert \\, Z) p(Z)$, где в аналитическом виде даны и $p(X \\, \\vert \\, Z)$, и $p(Z)$. При этом $p(X \\, \\vert \\, Z)$ называется правдоподобием выборки $X$ при условии $Z$, а $p(Z)$ является априорным распределением $Z$ (взятым на основании знаний о предметной области или полученным как апостериорное распределение $Z$ при условии какой-то ранее наблюдавшейся выборки $X_\\mathrm{old}$).\n",
"\n",
"Решение задачи, данной в байесовской постановке, описывается так:\n",
"* Апостериорное распределение находится по теореме Байеса:\n",
"$$p(Z \\, \\vert \\, X) = \\frac{p(X, Z)}{p(X)} = \\frac{p(X, Z)}{\\int p(X, Z) dZ},$$\n",
"* Точечная оценка $\\hat{Z}$ получается как мода или как математическое ожидание полученного апостериорного распределения:\n",
" - $\\hat{Z} = \\arg \\max_Z p(Z \\, \\vert \\, X)$,\n",
" - $\\hat{Z} = \\mathbb{E}_{p(Z \\, \\vert \\, X)}[Z]$.\n",
" \n",
"Основная сложность, связанная с байесовским статистическим выводом, как правило, заключается во взятии интеграла из знаменателя правой части формулы для апостериорного распределения.\n",
"\n",
"#### Способы решения задачи\n",
"\n",
"Существуют следующие способы найти апостериорное распределение $p(Z \\, \\vert \\, X)$:\n",
"* Если это возможно, явно взять интеграл из знаменателя в выражении для апостериорного распределения. Данный подход примением, например, тогда, когда $Z$ принимает дискретные значения и интеграл, на самом деле, является суммой, или тогда, когда $p(X \\, \\vert \\, Z)$ и $p(Z)$ — пара сопряжённых (conjugate) распределений;\n",
"* [Вариационный байесовский вывод](__home_url__/notes/Вариационный байесовский вывод и mean-field approximation). Для произвольного зафиксированного $X$ проводится приближённое вычисление апостериорного распределения $p(Z \\, \\vert \\, X)$ методом покоординатного подъёма, где на всех итерациях приближения ищутся в классе распределений, факторизующихся на непересекающиеся группы переменных;\n",
"* [Методы Монте-Карло по схеме марковской цепи](__home_url__/notes/Методы Монте-Карло над марковскими цепями), MCMC. Это численные методы взятия интеграла из знаменателя в выражении для апостериорного распределения, позволяющие сэмплировать значения $Z$ из распределения $p(Z \\, \\vert \\, X)$."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы",
"статистический_вывод"
]
},
"source": [
"## Вариационный байесовский вывод и mean-field approximation\n",
"\n",
"#### Описание и обоснование\n",
"\n",
"Пусть рассматривается задача [байесовского статистического вывода](__home_url__/notes/Байесовский статистический вывод), где по кортежу наблюдаемых переменных $X$ требуется найти апостериорное распределение кортежа скрытых переменных $p(Z \\, \\vert \\, X)$, имея явный вид вероятностной модели $p(X, Z) = p(X \\, \\vert \\, Z) p(Z)$, где обе плотности из правой части даны аналитически. Более детальное описание задачи и, в частности, то, как интерпретируются $X$ и $Z$, можно найти в заметке про байесовский статистический вывод.\n",
"\n",
"Вариационный байесовский вывод применяется в ситуации, когда апостериорное распределение $p(Z \\, \\vert \\, X) = \\frac{p(X, Z)}{p(X)} = \\frac{p(X, Z)}{\\int p(X, Z) dZ}$ не получается найти явно из-за того, что не вычисляется аналитически $p(X) = \\int p(X, Z) dZ$. Вместо искомого распределения $p(Z \\, \\vert \\, X)$ ищется его некоторое приближение $q(Z)$, тоже являющееся вероятностным распределением кортежа $Z$ (причём для разных $X$, вообще говоря, будут найдены разные $q(Z)$). Итак, цель заключается в том, чтобы найти $q(Z) \\approx p(Z \\, \\vert \\, X)$, где близость понимается в смысле минимизации дивергенции Кульбака-Лейблера между распределениями $q(Z)$ и $p(Z \\, \\vert \\, X)$:\n",
"$$D_{KL}(q(Z) \\, \\Vert \\, p(Z \\, \\vert \\, X)) \\to \\min_{q},$$\n",
"причём порядок аргументов дивергенции Кульбака-Лейблера такой, какой есть, потому что именно для такой дивергенции Кульбака-Лейблера удаётся найти решение сформулированной оптимизационной задачи.\n",
"\n",
"Для решения поставленной задачи сведём её к эквивалентной ей оптимизационной задаче.\n",
"\n",
"Каким бы ни было распределение $q(Z)$, всегда можно переписать логарифм правдоподобия наблюдаемой части выборки следующим образом:\n",
"$$\\log p(X) = \\log p(X) \\int q(Z) dZ = \\int q(Z) \\log p(X) dZ = \\int q(Z) \\log \\frac{p(X, Z)}{p(Z \\, \\vert \\, X)} dZ = \\int q(Z) \\log \\frac{p(X, Z)q(Z)}{p(Z \\, \\vert \\, X)q(Z)} dZ \\\\ = \\int q(Z) \\log \\frac{p(X, Z)}{q(Z)} dZ + \\int q(Z) \\log \\frac{q(Z)}{p(Z \\, \\vert \\, X)} dZ.$$\n",
"\n",
"Обозначим первое слагаемое последнего выражения за $\\mathcal{L}(q)$ и заметим, что второе слагаемое является той самой дивергенцией Кульбака-Лейблера $D_{KL}(q(Z) \\, \\Vert \\, p(Z \\, \\vert \\, X))$, которую необходимо минимизировать. Значит,\n",
"$$\\log p(X) = \\mathcal{L}(q) + D_{KL}(q(Z) \\, \\Vert \\, p(Z \\, \\vert \\, X)).$$\n",
"\n",
"Благодаря такому равенству можно взглянуть на поиск приближения $q(Z)$ для неизвестного распределения $p(Z \\, \\vert \\, X)$, которое, собственно говоря, неизвестно из-за невозможности вычислить $p(X)$, с двух точек зрения:\n",
"* Исходная оптимизационная задача. Минимизацией $D_{KL}(q(Z) \\, \\Vert \\, p(Z \\, \\vert \\, X))$ по $q$ можно получить приближение для $p(Z \\, \\vert \\, X)$, потому что дивергенция Кульбака-Лейблера неотрицательна и, более того, обнуляется тогда и только тогда, когда два распределения идентичны.\n",
"* Эквивалентная ей оптимизационная задача. Максимизацией $\\mathcal{L}(q)$ по $q$ можно получить приближение (а именно оценку снизу) для $p(X)$, потому что $\\mathcal{L}(q)$ не может быть больше $\\log p(X)$ в силу неотрицательности дивергенции Кульбака-Лейблера;\n",
"\n",
"Раз две выписанных оптимизационных задачи эквивалентны (имеют одинаковое решение), остановимся на второй:\n",
"$$\\mathcal{L}(q) \\to \\max_{q}.$$\n",
"Она выбрана, потому что удобнее для применения вариационного подхода. Функционал $\\mathcal{L}(q)$ называется нижней вариационной оценкой (evidence lower bound, ELBO), и он имеет также следующую интерпретацию:\n",
"$$\\mathcal{L}(q) = \\int q(Z) \\log(p(X \\, \\vert \\, Z)p(Z)) - \\log q(Z) \\, dZ = \\int q(Z) \\log p(X \\, \\vert \\, Z) - \\log \\frac{q(Z)}{p(Z)} \\, dZ = \\mathbb{E}_{Z \\sim q(Z)}[\\log p(X \\, \\vert \\, Z)] - D_{KL}(q(Z)\\, \\Vert \\, p(Z)),$$\n",
"что с точки зрения общей теории машинного обучения выглядит как сумма функционала качества и регуляризатора, штрафующего отклонения оценки $q(Z)$ апостериорного распределения $p(Z \\, \\vert \\, X)$ от априорного распределения $p(Z)$.\n",
"\n",
"Заметим, что до сих пор не вводилось никаких дополнительных предположений по сравнению с теми, что есть в общей постановке задачи байесовского статистического вывода. Вовлечение вариационного подхода опирается на то, что на распределение $q$ накладывается ограничение факторизуемости. Оно означает, что существует такая разбивка элементов кортежа $Z$ на $m$ непересекающихся множеств $s_i$, что\n",
"$$q(Z) = \\prod_{i=1}^m q_i(Z_{s_i}),$$\n",
"где $Z_{s_i}$ — кортеж, получающийся из кортежа $Z$ оставлением тех и только тех элементов, которые попадают в множество $s_i$, а $q_i(Z_i)$ — произвольное распределение кортежа $Z_{s_i}$. Стало быть, теперь ищется приближение распределения $p(Z \\, \\vert \\, X)$ в классе факторизуемых вероятностных распределений.\n",
"\n",
"Будем решать оптимизационную задачу итеративно методом покоординатного подъёма. Инициализируем как-то все $q_i$, а затем на каждом шаге будем по очереди обновлять каждое из них в соответствии со следующим правилом:\n",
"$$q_i(Z_{s_i}) := \\frac{\\exp(\\int \\log p(X, Z) \\prod_{j \\ne i} q_j(Z_{s_j}) dZ_{s_{-i}})}{\\int \\exp(\\int \\log p(X, Z) \\prod_{j \\ne i} q_j(Z_{s_j}) dZ_{s_{-i}}) dZ_{s_i}},$$\n",
"где $Z_{s_{-i}}$ обозначает кортеж, получающийся из $Z$ выкидыванием элементов, принадлежащих множеству $s_i$. Выписанная формула обновления $q_i(Z_{s_i})$ называется основной формулой вариационного байесовского вывода, а её получение здесь ради краткости опущено. Вывод же назван вариационным, потому что при оптимизации используется взятие вариации оптимизируемого функционала $\\mathcal{L}(q)$ вдоль направления, соответствующего $i$-му множеству скрытых переменных. На самом деле, существуют и другие подходы к вариационному байесовскому выводу, поэтому для большей конкретики описанный здесь подход называют mean-field approximation. Название пришло из физики магнитных полей, где есть раздел, именуемый теория среднего поля.\n",
"\n",
"Может показаться, что интеграл из знаменателя правой части основной формулы вариационного байесовского вывода ничуть не проще исходного интеграла из формулы для апостериорного распределения. Однако взять этот интеграл иногда удаётся. В следующем разделе будет описано, что и при каких условиях возможно.\n",
"\n",
"#### Применимость в зависимости от имеющихся свойств\n",
"\n",
"Чтобы описывать различные ситуации, введём ряд определений.\n",
"\n",
"Начнём с понятия условной сопряжённости, являющегося ослаблением понятия сопряжённости. Пусть распределение кортежа $Z$ факторизуется в разбивке элементов кортежа на $m$ непересекающихся множеств. Примем за $\\{Z_{s_i}\\}_{i=1}^m$ множество из $m$ кортежей, $i$-й из которых получается из $Z$ оставлением только соответствующих $i$-му множеству элементов, а через $Z_{s_{-i}}$ обозначим кортеж, получающийся из $Z$ выкидыванием всех элементов, соответствующих $i$-му множеству. Будем говорить, что присутствует условная сопряжённость, если одновременно выполнены $m$ условий, $i$-е из которых имеет вид:\n",
"* Для любого зафиксированного значения кортежа $Z_{s_{-i}}$ (то есть для любых зафиксированных значений всех кортежей $Z_{s_j}$, где берутся все $j$ от 1 до $m$ кроме $i$), распределения $p(X \\, \\vert \\, Z_{s_i}, Z_{s_{-i}})$ и $p(Z_{s_i} \\, \\vert \\, Z_{s_{-i}})$ сопряжены друг к другу.\n",
"\n",
"Теперь предположим, что кортеж $Z$ может быть разделён на два кортежа, обозначаемые далее как $T$ и $\\Theta$. Будем считать, что $p(Z) = p(T, \\Theta) = p(T)p(\\Theta)$, то есть имеется факторизация при разбивке элементов $Z$ на те, которые относятся к $T$, и на те, которые относятся к $\\Theta$. Для интерпретируемости далее можно считать, что кортеж $T$ составлен из строк матрицы скрытых признаков объектов, а $\\Theta$ — кортеж неизвестных параметров вероятностной модели.\n",
"\n",
"Определим $T$-сопряжённость как ослабление условной сопряжённости до того, что при факторизации $Z$ в разбивке на $T$ и $\\Theta$ условие выполнено только для $Z_{s_i} = T$, то есть как то, что пару сопряжённых распределений образуют $p(X \\, \\vert T, \\Theta)$ и $p(T \\, \\vert \\, \\Theta)$. Аналогично, определим $\\Theta$-сопряжённость как то, что сопряжены друг к другу распределения $p(X \\, \\vert \\Theta, T)$ и $p(\\Theta \\, \\vert \\, T)$.\n",
"\n",
"Далее, ослабим понятия $T$-сопряжённости и $\\Theta$-сопряжённости. Назовём условной $T$-сопряжённостью свойство, при котором распределение кортежа $T$ само как-то факторизуется и для факторизации распределения всего $Z$ на $\\Theta$ и на то, на что разбивается $T$, выполнены все условия из определения условной сопряжённости кроме того условия, где фиксируется всё, на что разбивается $T$ (то есть условия $\\Theta$-сопряжённости). Аналогично, условной $\\Theta$-сопряжённостью назовём свойство, при котором распределение кортежа $\\Theta$ как-то факторизуется и для факторизации распределения всего $Z$ на $T$ и на то, на что разбивается $\\Theta$, есть все условия из определения условной сопряжённости кроме того условия, где фиксируется всё, на что разбивается $\\Theta$ (то есть условия $T$-сопряжённости). \n",
"\n",
"В зависимости от присутствия того или иного свойства применимы разные методы статистического вывода. При этом с установлением того, какое свойство присутствует, есть одна тонкость: факторизуемости может и не быть, но всегда можно предполагать, будто она есть. После такого предположения будет искаться лишь приближение в классе соответствующим образом факторизованных распределений, и поэтому чем меньше искусственно разбивается кортеж $Z$, тем выше точность вывода (с другой стороны, тем выше и вычислительные затраты).\n",
"\n",
"Итак, соответствие между присутствующими свойствами и доступными методами вывода таково:\n",
"* сопряжённость (полная, т.е. $Z$-сопряжённость): полный байесовский вывод, без каких бы то ни было приближений оценивающий $p(Z \\, \\vert X) = p(T, \\Theta \\, \\vert \\, X)$;\n",
"* условная сопряжённость в разбивке $Z$ на $T$ и $\\Theta$: вариационный байесовский вывод, оценивающий $q(Z) = q(T)q(\\Theta) \\approx p(T, \\Theta \\, \\vert \\, X)$;\n",
"* $T$-сопряжённость: точный EM-алгоритм, являющийся частным случаем вариационного байесовского вывода, где дополнительно предполагается, что $q(\\Theta) = \\delta(\\Theta - \\Theta_0)$ и в этой записи $\\delta$ — дельта-функция Дирака, а $\\Theta_0$ — константный кортеж, который требуется найти; точный EM-алгоритм оценивает $q(Z) = q(T)\\delta(\\Theta - \\Theta_0) \\approx p(T, \\Theta \\, \\vert \\, X)$;\n",
"* $\\Theta$-сопряжённость: точный ME-алгоритм, являющийся частным случаем вариационного байесовского вывода, где дополнительно предполагается, что $q(T) = \\delta(T - T_0)$, а $T_0$ — константный кортеж, который требуется найти; точный ME-алгоритм оценивает $q(Z) = \\delta(T - T_0)q(\\Theta) \\approx p(T, \\Theta \\, \\vert \\, X)$;\n",
"* условная $T$-сопряжённость: вариационный EM-алгоритм, оценивающий $q(Z) = \\left(\\prod_{i=1}^m q_i(T_i)\\right) \\delta(\\Theta - \\Theta_0) \\approx p(T, \\Theta \\, \\vert \\, X)$;\n",
"* условная $\\Theta$-сопряжённость: вариационный ME-алгоритм, оценивающий $q(Z) = \\delta(T - T_0) \\left(\\prod_{i=1}^m q_i(\\Theta_i)\\right) \\approx p(T, \\Theta \\, \\vert \\, X)$;\n",
"* никаких свойств нет: жёсткий EM-алгоритм, оценивающий $q(Z) = \\delta(T - T_0) \\delta(\\Theta - \\Theta_0) \\approx p(T, \\Theta \\, \\vert \\, X)$, то есть дающий только точечные оценки.\n",
"\n",
"Получившийся список можно прокомментировать так. Чем более слабое свойство имеется, тем более «простым» распределением приходится полагать то распределение, которое ищется в задаче статистического вывода. «Простота» здесь понимается как факторизуемость и атомарность.\n",
"\n",
"Напоследок разберём какой-нибудь пример, показывающий, почему указанных свойств достаточно для соответствующих методов. Скажем, остановимся на достаточности $T$-сопряжённости для точного EM-алгоритма. В точном EM-алгоритме сложное место находится на E-шаге, где ищется распределение $$q^{(n+1)}(Z) = p(Z \\, \\vert \\, X, \\Theta^{(n)}) = \\frac{p(X, Z \\, \\vert \\, \\Theta^{(n)})}{\\int p(X, Z \\, \\vert \\, \\Theta^{(n)}) dZ} = \\frac{p(X \\, \\vert \\, Z, \\Theta^{(n)})p(Z \\, \\vert \\, \\Theta^{(n)})}{\\int p(X \\, \\vert \\, Z, \\Theta^{(n)})p(Z \\, \\vert \\, \\Theta^{(n)}) dZ} = \\frac{p(X \\, \\vert \\, T, \\Theta^{(n)})p(T \\, \\vert \\, \\Theta^{(n)})\\delta(\\Theta - \\Theta^{(n)})}{\\int p(X \\, \\vert \\, T, \\Theta^{(n)})p(T \\, \\vert \\, \\Theta^{(n)})\\delta(\\Theta - \\Theta^{(n)}) dT d\\Theta} = \\frac{p(X \\, \\vert \\, T, \\Theta^{(n)})p(T \\, \\vert \\, \\Theta^{(n)})}{\\int p(X \\, \\vert \\, T, \\Theta^{(n)})p(T \\, \\vert \\, \\Theta^{(n)}) dT} = \\frac{p(X, T \\, \\vert \\, \\Theta^{(n)})}{\\int p(X, T \\, \\vert \\, \\Theta^{(n)}) dT},$$\n",
"и тут интеграл в полученном в конце выражении берётся явно в силу $T$-сопряжённости."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы",
"MCMC"
]
},
"source": [
"## Методы Монте-Карло над марковскими цепями\n",
"\n",
"#### Введение\n",
"\n",
"Пусть с точностью до константы известна в аналитическом виде плотность $p(x)$ распределения случайного вектора $x$. Требуется научиться сэмплировать выборки из этого распределения, то есть получать объекты $x_i \\sim p(x)$, или хотя бы научиться приближённо интегрировать по этому распределению, то есть оценивать интегралы $\\int f(x)p(x) dx$ для любой достаточно регулярной функции $f$. Решение второй задачи получается из решения первой задачи, ведь выписанный интеграл аппроксимируется средним $f(x)$ по просэмплированной выборке достаточного размера.\n",
"\n",
"Умение решать так поставленную задачу пригождается, в частности, при проведении [байесовского статистического вывода](__home_url__/notes/Байесовский статистический вывод). Апостериорное распределение $p(Z \\, \\vert \\, X)$ может быть получено в аналитическом виде с точностью до нормировочной константы, и затем при помощи сэмплирования из него можно будет по-байесовски ансамблировать предсказания, выдаваемые моделями с различными значениями неизвестных параметров.\n",
"\n",
"В случае, когда $p(x)$ является распределением над низкоразмерным вектором, существуют алгоритмы, позволяющие сэмплировать независимые объекты в точности из распределения $p(x)$. К таким алгоритмам относятся сэмплирование с отклонением (rejection sampling) и сэмплирование с перевзвешиванием (importance sampling). Однако оба упомянутых алгоритма подвержены проклятию размерности, так что для данных высокой размерности нужны другие методы. Далее речь пойдёт о методах, которые, с одной стороны, более устойчивы к проклятью размерности, но, с другой стороны, сэмплируют так, что результаты соседних сэмплирований не являются независимыми (этот недостаток можно смягчить, оставляя лишь результаты каждого $m$-го испытания, где $m$ достаточно велико).\n",
"\n",
"#### Схема Метрополиса-Хастингса\n",
"\n",
"Введём следующие сущности:\n",
"* $\\overline{p}(x)$ — функция, с точностью до константы равная функции плотности $p(x)$ интересующего вероятностного распределения;\n",
"* $q(x \\, \\vert \\, y)$, где $y \\in \\mathrm{Domain}(x)$ — некоторая плотность вероятностного распределения, называемого предложным распределением и обладающего такими свойствами:\n",
" - $\\forall y$ есть возможность сэмплировать из распределения $q(x \\, \\vert \\, y)$,\n",
" - $\\forall x, y$, таких что $p(x) > 0$ и $p(y) > 0$, также и $q(x \\, \\vert \\, y) > 0$.\n",
" \n",
"В качестве предложного распределения может быть использовано, например, $q(x \\, \\vert \\, y) = \\mathcal{N}_{y, \\sigma^2I}(x)$. Такое распределение симметрично, $q(x \\, \\vert \\, y) = q(y \\, \\vert \\, x)$, и поэтому с ним может работать даже схема Метрополиса, то есть предшественница схемы Метрополиса-Хастингса, которую Хастингс обобщил на случай несимметричных предложных распределений.\n",
"\n",
"Идея заключается в том, что хочется построить эргодическую марковскую цепь, инвариантным распределением которой будет $p(x)$. Если такая цепь будет построена, начиная с какого-то момента сэмплировать из $p(x)$ можно будет за счёт порождения следующего элемента этой цепи.\n",
"\n",
"Сначала опишем конструкцию схемы Метрополиса-Хастингса, а потом докажем, что она работает так, как заявлено.\n",
"\n",
"Итак, пусть текущий элемент цепи равен $x_n$. Следующий за ним элемент цепи $x_{n+1}$ порождается по представленному ниже алгоритму:\n",
"* Сэмплируется $x \\sim q(x \\, \\vert x_n)$;\n",
"* Вычисляется отношение Хастингса $H(x, x_n) = \\min\\left(1, \\frac{\\overline{p}(x)q(x_n \\, \\vert \\, x)}{\\overline{p}(x_n)q(x \\, \\vert \\, x_n)}\\right) = \\min\\left(1, \\frac{p(x)q(x_n \\, \\vert \\, x)}{p(x_n)q(x \\, \\vert \\, x_n)}\\right)$;\n",
"* С вероятностью, равной $H(x, x_n)$, $x_{n+1} := x$, а иначе $x_{n+1} := x_n$.\n",
"\n",
"Начало цепи $x_1$ может быть выбрано произвольно, но от его выбора может зависеть, как долго придётся ждать, пока марковская цепь сойдётся к инвариантному распределению.\n",
"\n",
"Теперь перейдём к доказательству работоспособности конструкции.\n",
"\n",
"Существует теорема, гласящая, что любая однородная марковская цепь с переходной вероятностью $p(x \\, \\vert \\, y)$, такой что $p(x \\, \\vert \\, y) > 0$ для всех $x$ и $y$, является эргодической, то есть имеет ровно одно инвариантное распределение, к которому непременно сходится. В случае со схемой Метрополиса-Хастингса при $x_{n+1} \\ne x_n$ переходная вероятность равна $q(x_{n+1} \\, \\vert \\, x_{n})H(x_{n+1}, x_n) > 0$, и при $x_{n+1} = x_n$ она тоже больше нуля, так что эта теорема применима.\n",
"\n",
"Осталось показать, что инвариантное распределение имеет плотность $p(x)$. Чтобы сделать это, прибегнем к тривиальному факту, гласящему, что для марковской цепи с переходной вероятностью $r(x \\, \\vert \\, y)$ из того, что $r(y \\, \\vert \\, x)p(x) = r(x \\, \\vert y)p(y)$, следует, что $p(x)$ является инвариантным распределением. Предпосылку сделанного утверждения называют уравнением детального баланса, это уравнение говорит, что, грубо говоря, сколько вероятности переходит из $x$ в $y$ при переходе к следующему элементу марковской цепи, столько же переходит и из $y$ в $x$. Но это и означает, что $p(x)$ задаёт инвариантное распределение. Итак, убедимся, что для схемы Метрополиса-Хастингса уравнение детального баланса для $p(x)$ соблюдается (строго говоря, выкладка, приведённая ниже, применима только к случаю $x_{n+1} \\ne x_n$, но случай $x_{n+1} = x_n$ тривиален):\n",
"$$r(x_{n+1} \\, \\vert \\, x_n)p(x_n) = q(x_{n+1} \\, \\vert \\, x_{n})H(x_{n+1}, x_n)p(x_n) = p(x_n)q(x_{n+1} \\, \\vert \\, x_{n}) \\min\\left(1, \\frac{p(x_{n+1})q(x_n \\, \\vert \\, x_{n+1})}{p(x_n)q(x_{n+1} \\, \\vert \\, x_n)}\\right) = \\min\\left(p(x_n)q(x_{n+1} \\, \\vert \\, x_{n}), p(x_{n+1})q(x_n \\, \\vert \\, x_{n+1})\\right) = p(x_{n+1})q(x_n \\, \\vert \\, x_{n+1})\\min\\left(1, \\frac{p(x_n)q(x_{n+1} \\, \\vert \\, x_n)}{p(x_{n+1})q(x_n \\, \\vert \\, x_{n+1})}\\right) = r(x_n \\, \\vert \\, x_{n+1})p(x_{n+1}).$$\n",
"Тем самым выше было показано, что для построенной марковской цепи среди инвариантных распределений есть $p(x)$, но поскольку инвариантное распределение ровно одно, только $p(x)$ и может быть им.\n",
"\n",
"Неформально говоря, $H(x, x_n)$ для того и вводится, чтобы добиться соблюдения уравнения детального баланса. Также можно интерпретировать $H(x, x_n)$ как то, что побуждает марковскую цепь реже уходить из какой-либо точки, если в новой точке $p(x)$ ниже.\n",
"\n",
"#### Схема Гиббса\n",
"\n",
"Схема Гиббса является частным случаем схемы Метрополиса-Хастингса, имеющим то свойство, что всегда $H(x, x_n) = 1$. Достигается это за счёт того, что теперь шаг от $x_n$ к $x_{n+1}$ делается поэтапно, а именно для каждой координаты по порядку начиная с первой сначала проводится её обновление, затем происходит обновление следующей координаты, и так далее вплоть до последней координаты включительно.\n",
"\n",
"Если обозначить за $x_n^i$ $i$-ю компоненту вектора $x_n$, то тогда формула для её обновления в схеме Гиббса имеет вид:\n",
"$$x_{n+1}^k \\sim \\overline{p}(x^k \\, \\vert \\, x^1 = x_{n+1}^1, \\, \\ldots, \\, x^{k-1} = x_{n+1}^{k-1}, \\, x^{k+1} = x_n^{k+1}, \\, \\ldots, \\, x^d = x_n^d),$$\n",
"где $d$ — размерность векторов $x_n$ и $x_{n+1}$. Знак $\\sim$ в формуле обозначает, что $x_{n+1}^k$ сэмплируется из распределения, стоящего справа от этого знака. Сэмплировать из такого распределения можно, потому что оно является распределением одномерной случайной величины, известным с точностью до константы. Обоснуем это утверждение. Во-первых, с точностью до константы это распределение известно, потому что оно является ограничением известного с точностью до константы распределения $\\overline{p}$. Во-вторых, для известных с точностью до константы распределений одномерных случайных величин можно применять техники, такие как сэмплирование с перевзвешиванием, потому что проклятие размерности нестрашно.\n",
"\n",
"При предложном распределении $q(x^k \\, \\vert \\, y) = \\overline{p}(x^k \\, \\vert \\, y^1, ..., y^{k-1}, y^{k+1}, ..., y^d)$ выкладка показывает, что $H(x, x_n) = 1$. Если выписанное предложное распределение развернуть от распределения одной $k$-й компоненты до распределения всего $d$-мерного вектора, то и получится, что схема Гиббса является схемой Метрополиса-Хастингса, в которой выбрано развёрнутое предложное распределение."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"MCMC",
"байесовские_методы"
]
},
"source": [
"## Масштабирование марковских цепей Монте-Карло и схема Ланжевена\n",
"\n",
"#### Введение\n",
"\n",
"Пусть решается задача [байесовского статистического вывода](__home_url__/notes/Байесовский статистический вывод) и требуется оценить апостериорное распределение\n",
"$$p(Z \\, \\vert \\, X) = \\frac{p(X \\, \\vert \\, Z)p(Z)}{\\int p(X \\, \\vert \\, Z)p(Z) dZ}.$$\n",
"Числитель по предположению задачи байесовского статистического вывода известен в аналитическом виде, а вот знаменатель является константой, вычислить которую часто бывает невозможно из-за того, что интеграл не берётся. Приходится прибегать к численным методам, одним из которых являются [марковские цепи Монте-Карло](__home_url__/notes/Методы Монте-Карло над марковскими цепями). Они позволяют сэмплировать $Z$ из $p(Z \\, \\vert \\, X)$.\n",
"\n",
"Недостатком классических схем для марковских цепей Монте-Карло, то есть схемы Метрополиса-Хастингса и схемы Гиббса, является то, что они работают слишком медленно при большом количестве объектов в выборке (при большом количестве строк в матрице $X$).\n",
"\n",
"Так, в схеме Метрополиса-Хастингса отношение Хастингса для некоторого предложного распределения $q(Z \\, \\vert \\, Z_n)$ и просэмплированного из него кандидата $Z$ имеет вид:\n",
"$$H(Z, Z_n) = \\min\\left(1, \\frac{p(X \\, \\vert \\, Z)p(Z) \\, q(Z_n \\, \\vert \\, Z)}{p(X \\, \\vert \\, Z_n)p(Z_n) \\, q(Z \\, \\vert \\, Z_n)}\\right).$$\n",
"Тут плохо масштабируется вычисление $p(X \\, \\vert \\, Z)$ и $p(X \\, \\vert \\, Z_n)$. Скажем, если все объекты независимы, придётся перемножить столько слагаемых, сколько строк в $X$. Поскольку подобное перемножение делается при сэмплировании каждого нового значения $Z_{n+1}$, это существенно замедляет схему.\n",
"\n",
"В схеме же Гиббса проблема с масштабированием возникает из-за того, что каждая компонента вектора $Z_{n+1}$ сэмплируется из распределения, для вычисления которого с точностью до константы тоже требуется обработать всю матрицу $X$.\n",
"\n",
"Существуют модификации марковских цепей Монте-Карло с улучшенной масштабируемостью на большие данные, отталкивающиеся от следующей идеи. Можно заметить, что схемы Метрополиса-Хастингса и Гиббса пытаются найти $Z$, для которых высока $p(Z \\, \\vert \\, X)$, используя для этого только значения $p(Z \\, \\vert \\, X)$. Однако что, если также использовать градиент $p(Z \\, \\vert \\, X)$ или, точнее, градиент $\\log p(Z \\, \\vert \\, X)$, в случае независимых объектов распадающийся на сумму градиентов? Если так поступить, вместо градиента можно будет использовать стохастический градиент, оцениваемый только по текущему батчу, что и снимет проблему большого размера $X$.\n",
"\n",
"#### Схема Ланжевена\n",
"\n",
"Схема Ланжевена опирается на концепцию из физики, называемую динамика Ланжевена (Langevin dynamics). На базе такой динамики может быть постоена схема Метрополиса-Хастингса с определённым предложным распределением. Сама по себе эта схема не является масштабируемой, но её можно переделать в масштабируемую по аналогии с тем, как из классического градиентного подъёма получается стохастический градиентный подъём.\n",
"\n",
"Итак, сначала опишем схему Метрополиса-Хастингса на базе динамики Ланжевена. Чтобы её описать, достаточно задать предложное распределение. Зададим его так:\n",
"$$q(Z \\, \\vert \\, Z_n) = \\mathcal{N}_{\\mu, \\varepsilon}(Z),$$\n",
"где $\\mu = Z_n + \\frac{\\varepsilon}{2}\\nabla_Z\\left(\\log(p(X \\, \\vert \\, Z_n)p(Z_n))\\right)$, $\\varepsilon \\in \\mathbb{R}$, а $\\nabla_Z$ обозначает градиент по $Z$ (в данном случае, вычисленный в точке $Z_n$). Подчеркнём, что $\\varepsilon$ фигурирует и как дисперсия предложного распределения, и как множитель перед градиентом в выражении для среднего. Также отметим, что выписанное предложное распределение не является симметричным, и поэтому с ним работает именно схема Метрополиса-Хастингса, а не схема Метрополиса.\n",
"\n",
"Теперь представим, что почему-то хочется проигнорировать коррекцию Метрополиса-Хастингса и, не вычисляя отношение Хастингса, всегда принимать просэмплированные из предложного распределения объекты. Тогда получится следующая схема (марковская цепь):\n",
"$$Z_{n+1} = Z_n + \\frac{\\varepsilon}{2}\\nabla_Z\\left(\\log(p(X \\, \\vert \\, Z_n)p(Z_n))\\right) + \\eta_{n},$$\n",
"где $\\eta_{n} \\sim \\mathcal{N}_{0, \\varepsilon}(\\eta_{n})$. Если бы не последнее слагаемое, это была бы формула поиска моды распределения $p(Z \\, \\vert \\, X)$ методом покоординатного подъёма, ведь\n",
"$$\\nabla_Z\\left(\\log p(Z_n \\, \\vert \\, X)\\right) = \\nabla_Z\\left(\\log(p(X \\, \\vert \\, Z_n)p(Z_n)) - \\log(\\int p(X \\, \\vert \\, Z)p(Z) dZ)\\right) = \\nabla_Z\\left(\\log(p(X \\, \\vert \\, Z_n)p(Z_n))\\right).$$\n",
"Однако благодаря добавлению шума, вносимого слагаемым $\\eta_{n}$, схема описывает некоторое блуждание, как и в случае схемы Метрополиса-Хастингса склонное задерживаться в регионах с высокой плотностью $p(Z \\, \\vert \\, X)$.\n",
"\n",
"Через стохастические дифференциальные уравнения можно доказать, что если константу $\\varepsilon$ заменить на элемент последовательности $\\varepsilon_n$, такой что $\\varepsilon_n \\rightarrow 0$, то тогда описанная марковская цепь тоже будет порождать $Z_{n}$ из распределения $p(Z \\, \\vert \\, X)$. Достигается это за счёт того, что множитель перед градиентом и дисперсия шума соотносятся друг с другом (а именно первое в два раза меньше второго).\n",
"\n",
"Избавление от необходимости вычислять отношение Хастингса ускоряет вычисления — правда, ценой того, что из-за стремящихся к нулю $\\varepsilon_n$ с какого-то момента соседние просэмплированные значения практически не будут отличаться друг от друга (т.е. будут сильно скоррелированы), из-за чего потребуется прореживание, становящееся всё более сильным.\n",
"\n",
"Если объекты выборки независимы друг от друга, то можно получить второй существенный выигрыш. Разобъём выборку на батчи размера $b$, $n$-й из которых имеет вид $x_{n_1}$, ..., $x_{n_b}$. Тогда рассмотрим такую схему:\n",
"$$Z_{n+1} = Z_n + \\frac{\\varepsilon_n}{2}\\left(\\frac{N}{b}\\sum_{j=1}^b \\nabla_Z\\log(p(x_{n_j} \\, \\vert \\, Z_n)p(Z_n))\\right) + \\eta_{n},$$\n",
"где $N$ — количество строк в матрице $X$, а $\\eta_{n} \\sim \\mathcal{N}_{0, \\varepsilon_n}(\\eta_{n})$. По построению новая схема на каждом шаге использует лишь батч ограниченного размера, а не всю выборку.\n",
"\n",
"Чтобы новая схема сэмплировала $Z$ из распределения $p(Z \\, \\vert \\, X)$, достаточно наложить такие ограничения на последовательность $\\varepsilon_n$: $\\sum_{n=1}^{+\\infty} \\varepsilon_n = +\\infty$ и $\\sum_{n=1}^{+\\infty} \\varepsilon_n^2 < +\\infty$. В качестве примера подходящей последовательности подойдёт, например, $\\varepsilon_n = \\frac{1}{n^d}$, где $0.5 < d \\le 1$. Доказательство того, что наложенных ограничений достаточно, опирается на то, что при вычислении $Z_{n+1}$ есть два источника случайности: выбор объектов в батч и добавление $\\eta_{n}$. В силу наложенных ограничений дисперсия, вносимая первым из них, со временем станет пренебрежимо малой в сравнении с дисперсией, вносимой вторым из них. Грубо говоря, если оставить только второй источник случайности, то получится то же самое, что было в предыдущей схеме, то есть в схеме без разбивки на батчи. А у такой схемы, как уже упоминалось, работоспособность доказывается через стохастические дифференциальные уравнения. \n",
"\n",
"#### Альтернативы\n",
"\n",
"Марковские цепи Монте-Карло можно строить и на динамике Гамильтона. Такая конструкция имеет то преимущество, что в случае мультимодального распределения $p(Z \\, \\vert \\, X)$ переключение между модами, возле которых сэмплируются значения $Z$, происходит чаще."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"графические_модели"
]
},
"source": [
"## Вероятностные графические модели\n",
"\n",
"#### Введение\n",
"\n",
"Рассмотрим совместное распределение $n$ случайных величин $\\{x_i\\}_{i=1}^n$. Многкратно применяя определение условной вероятности, это распределение можно записать в таком виде:\n",
"$$p(\\{x_i\\}_{i=1}^n) = \\prod_{j=1}^n p(x_j \\vert \\{x_i\\}_{i=1}^{j-1}),$$\n",
"где при $j=1$ условная вероятность $x_1$ при условии пустого множества $\\{x_i\\}_{i=1}^{0}$ понимается просто как $p(x_1)$. Тем самым исходное вероятностное распределение, каким бы оно ни было, представлено в виде произведения условных вероятностей.\n",
"\n",
"Если $n$ хоть сколько-нибудь велико, оценивать $p(\\{x_i\\}_{i=1}^n)$ затруднительно как с вычислительной точки зрения, так и со статистической точки зрения («проклятие размерности»). Введя ряд дополнительных предположений о $p(\\{x_i\\}_{i=1}^n)$, его можно представить в виде произведения вероятностей, которое имело бы более удобный вид, чем выписанный выше.\n",
"\n",
"Вероятностная графическая модель некого распределения сопоставляет этому распределению граф с $x_i$ в качестве вершин. При помощи данного графа можно наглядно определить, какие множители участвуют в разложении $p(\\{x_i\\}_{i=1}^n)$. При этом выбор графа в некоторой степени условен. Одному и тому же распределению могут быть сопоставлены несколько разных графов, и одному и тому же графу могут соответствовать разные совместные распределения.\n",
"\n",
"#### Направленные графические модели\n",
"\n",
"Направленная графическая модель сопоставляет вероятностному распределению ориентированный граф без циклов (циклы здесь подразумеваются с учётом направления). Для ориентированного графа без циклов верно, что его вершины можно пронумеровать так, чтобы каждое ребро исходило из вершины с меньшим номеров, чем у той вершины, куда оно входит. Будем считать, что случайные величины $\\{x_i\\}_{i=1}^n$ проиндексированы подобным образом. Тогда общий случай, разобранный в предыдущем разделе, соответствует графу, где из вершины, соответствующей $x_i$, исходят рёбра в вершины, соответствующие всем $x_j$, где $j > i$.\n",
"\n",
"Упрощение разложения достигается тогда, когда в графе отсутствует часть рёбер по сравнению с графом для общего случая. Вероятностная графическая модель для такого графа предполагает, что совместное распределение имеет следующий вид:\n",
"$$p(\\{x_i\\}_{i=1}^n) = \\prod_{j=1}^n p(x_j \\vert P_j),$$\n",
"где за $P_j$ обозначено множество случайных величин, таких что из соответствующих им вершин исходят рёбра в $j$-ю вершину.\n",
"\n",
"#### Ненаправленные графические модели\n",
"\n",
"Ненаправленная графическая модель сопоставляет вероятностному распределению обычный (то есть не ориентированный) граф, который может содержать циклы. В этом графе вершины, соответствующие $x_i$ и $x_j$, соединены ребром тогда и только тогда, когда $x_i$ и $x_j$ не являются условно независимыми при условии всех остальных случайных величин. Таким образом, из полносвязного графа, соответствующего общему случаю, когда про совместное распределение ничего неизвестно, удалением рёбер можно получить граф, по которому будет построено более удобное разложение на множители.\n",
"\n",
"В ненаправленных графических моделях разложение ищется в следующем виде:\n",
"$$p(\\{x_i\\}_{i=1}^n) = \\frac{1}{z} \\prod_{j=1}^{l} \\phi_j(S_j),$$\n",
"где $S_j$ — множества, составленные из исходных случайных величин $\\{x_i\\}_{i=1}^n$, такие, что всего их $l$ штук, а ещё верно, что $\\bigcup_{j=1}^l S_j = \\{x_i\\}_{i=1}^n$, $\\phi_j$ — произвольные функции со значениями в $\\mathbb{R}_{\\ge 0}$, а $z$ — нормировочная константа, благодаря которой полная вероятность равна 1. Нелишне подчеркнуть, что множества $S_j$ могут пересекаться друг с другом.\n",
"\n",
"Чтобы понять, в каком виде можно брать $S_j$, нужно ввести несколько понятий из теории графов. Кликой называется множество вершин графа, такое что любые две вершины из него соединены ребром. Максимальной кликой называется такая клика, что она не содержится в какой-либо другой клике.\n",
"\n",
"Для графа, сопоставленного совместному распределению случайных величин $\\{x_i\\}_{i=1}^n$, обозначим количество его максимальных клик за $m$, а сами максимальные клики — за $C_j$, $1 \\le j \\le m$. Утверждается, что без потери общности разложение можно искать в таком виде:\n",
"$$p(\\{x_i\\}_{i=1}^n) = \\frac{1}{z} \\prod_{j=1}^{m} \\phi_j(C_j).$$\n",
"\n",
"Чтобы доказать это, достаточно проверить, что:\n",
"* нет смысла рассматривать множество $S_j$, где какие-либо две случайных величины условно независимы при условии всех остальных величин,\n",
"* нет смысла рассматривать немаксимальные клики.\n",
"\n",
"Проверим первое утверждение. Пусть $x_i \\in S_j$ и $x_k \\in S_j$, но между $i$-й и $k$-й вершинами нет ребра. Обозначим за $X_{-\\{i, k\\}}$ множество оставшихся случайных величин. Тогда:\n",
"$$p(x_i, x_k \\vert X_{-\\{i, k\\}}) = p(x_i \\vert X_{-\\{i, k\\}}) p(x_k \\vert X_{-\\{i, k\\}}) \\propto \\phi_j(\\{x_i\\} \\cup \\{x_k\\} \\cup (S_j \\cap X_{-\\{i, k\\}})),$$\n",
"где правая часть является $\\phi_j(S_j)$, переписанной так, чтобы выделить, что по $X_{-\\{i, k\\}}$ происходит обусловливание. Полученная пропорциональность говорит о том, что $\\phi_j(S_j)$ можно разбить на множители аналогичного вида, каждый из которых не будет включать одновременно и $x_i$, и $x_k$. Таким образом, факторизации, где какие-либо две случайных величины условно независимы при условии всех остальных величин, допускают дальнейшую факторизацию.\n",
"\n",
"Что касается второго утверждения, то любой множитель, соответствующий немаксимальной клике, избыточен, если в разложении уже есть множитель, соответствующий содержащей её максимальной клике. Это так, потому что произведение функции нескольких аргументов и функции подмножества этих исходных аргументов тоже является функцией исходных аргументов. А для каждой максимальной клики должен быть множитель, потому что только он отражает все зависимости между случайными величинами, образующими эту клику."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"графические_модели",
"последовательности",
"EM"
]
},
"source": [
"## Скрытые марковские модели\n",
"\n",
"#### Введение\n",
"\n",
"Скрытая марковская модель с $m$ скрытыми состояниями и $n$ наблюдаемыми состояниями — это вероятностная графическая модель, описываемая совокупностью сущностей из следующего списка:\n",
"* Вероятностное распределение $\\pi$ на множестве $\\{1, ..., m\\}$, представляющее собой распределение скрытого состояния в начальный момент времени;\n",
"* Матрица переходов между скрытыми состояниями $A$, её размер $m \\times m$, а на пересечении $i$-й строки и $j$-го столбца стоит вероятность из скрытого состояния $i$ перейти в скрытое состояние $j$; как следствие сумма элементов любой строки матрицы $A$ равна 1;\n",
"* Матрица эмиссии наблюдаемых значений $B$, её размер $m \\times n$, а на пересечении $i$-й строки и $j$-го столбца стоит условная вероятность того, что текущее наблюдаемое состояние равно $j$ при условии, что текущее скрытое состояние есть $i$; сумма элементов любой строки матрицы $B$ тоже равна 1.\n",
"\n",
"Для переноса на случай, когда наблюдаемые значения непрерывны, можно считать, что вместо матрицы $B$ имеются $m$ вероятностных распределений.\n",
"\n",
"В контексте скрытых марковских моделей существуют следующие задачи:\n",
"* Зная $\\pi$, $A$ и $B$, вычислить вероятность данной последовательности наблюдаемых значений $(x_t)_{t=1}^T$, где $T$ — длина последовательности;\n",
"* Зная $\\pi$, $A$ и $B$, по имеющейся последовательности наблюдаемых значений $(x_t)_{t=1}^T$ восстановить последовательность скрытых значений $(h_t)_{t=1}^T$;\n",
"* Имея множество последовательностей наблюдаемых значений $\\left\\{\\left. (x_t)_{t=1}^{T_k} \\; \\right\\vert \\; k \\in \\{1, ..., N\\}\\right\\}$, где $N$ — количество обучающих последовательностей, оценить $\\pi$, $A$ и $B$, то есть обучить скрытую марковскую модель.\n",
"\n",
"#### Задача вычисления вероятности наблюдаемой последовательности\n",
"\n",
"Эта задача является задачей вероятностного (прямого) вывода, а не задачей статистического (обратного) вывода. На самом деле, подсчитать искомую вероятность можно просто в лоб по определению скрытой марковской модели:\n",
"$$p\\!\\left((x_t)_{t=1}^T\\right) = \\sum_{(h_t)_{t=1}^T \\in \\{1, ..., m\\}^T} \\left(\\pi(h_1) B_{h_1x_1} \\prod_{t=2}^T A_{h_{t-1}h_t} B_{h_tx_t}\\right),$$\n",
"однако такой подход имеет экспоненциальную от длины наблюдаемой последовательности $T$ сложность, потому что суммирование проводится по $m^T$ возможным вариантам последовательностей скрытых состояний.\n",
"\n",
"Для того чтобы задача была вычислительно подъёмной, используется алгоритм прохода вперёд (forward algorithm), восходящий к динамическому программированию.\n",
"\n",
"Введём индексируемые позицией в последовательности $t$ и скрытым состоянием $i$ переменные\n",
"$$\\alpha_{ti} = p\\!\\left((x_s)_{s=1}^t, h_t = i\\right).$$\n",
"Иными словами, $\\alpha_{ti}$ есть совместная вероятность появления подпоследовательности $(x_s)_{s=1}^t$ на первых $t$ наблюдаемых позициях и события, заключающегося в том, что в $t$-й позиции скрытое состояние равнялось $i$.\n",
"\n",
"Теперь можно описать алгоритм прохода вперёд:\n",
"1. Для всех $i \\in \\{1, ..., m\\}$ вычислить начальные значения:\n",
"$$\\alpha_{1i} = \\pi(i)B_{ix_1}.$$\n",
"2. Для всех $(t, i) \\in \\{2, ..., T\\} \\times \\{1, ..., m\\} $ в порядке неубывания $t$ провести индуктивный шаг вперёд:\n",
"$$\\alpha_{ti} = \\sum_{j=1}^m \\alpha_{(t-1)j}A_{ji}B_{ix_t}.$$\n",
"3. Вычислить искомую вероятность:\n",
"$$p\\!\\left((x_t)_{t=1}^T\\right) = \\sum_{i=1}^m \\alpha_{Ti}.$$\n",
"\n",
"Поскольку в алгоритме прохода вперёд возникают произведения длины $T$, все сомножители которых являются вероятностями, при хоть сколько-нибудь больших $T$ возникает вычислительная проблема округления слишком малых чисел до нуля (underflow). Для решения этой проблемы вводится нормализация, при которой вместо $\\alpha_{ti}$ подсчитываются отмасштабированные переменные $\\hat{\\alpha}_{ti}$, задаваемые следующими рекуррентными соотношениями:\n",
"$$\\hat{\\alpha}_{1i} = \\frac{\\pi(i)B_{ix_1}}{\\sum_{j=1}^m \\pi(j)B_{jx_1}},$$\n",
"$$\\hat{\\alpha}_{ti} = \\frac{\\sum_{j=1}^m \\hat{\\alpha}_{(t-1)j}A_{ji}B_{ix_t}}{\\sum_{k=1}^m \\left( \\sum_{j=1}^m \\hat{\\alpha}_{(t-1)j}A_{jk}B_{kx_t}\\right)}.$$\n",
"При фиксированном $t$ знаменатель из правой части для всех $\\hat{\\alpha}_{ti}$ один и тот же. Обозначим его за $c_t$. Тогда искомая вероятность равна $\\prod_{t=1}^T c_t$, потому что\n",
"$$\\prod_{t=1}^T c_t = \\prod_{t=1}^T c_t \\cdot \\sum_{i=1}^m \\hat{\\alpha}_{Ti} = \\prod_{t=1}^{T-1} c_t \\cdot \\sum_{i=1}^m \\left(\\sum_{j=1}^m \\hat{\\alpha}_{(T-1)j}A_{ji}B_{ix_T}\\right) = ... = \\sum_{i=1}^m \\alpha_{Ti} = p\\!\\left((x_t)_{t=1}^T\\right).$$\n",
"\n",
"#### Задача восстановления последовательности скрытых состояний\n",
"\n",
"Эта задача решается алгоритмом Витерби, который, как и алгоритм прохода вперёд, тоже восходит к динамическому программированию.\n",
"\n",
"Введём индексируемые позицией в последовательности $t$ и скрытым состоянием $i$ переменные\n",
"$$\\psi_{ti} = \\max_{(h_s)_{s=1}^t \\in \\{1, ..., m\\}^t}p\\!\\left((x_s)_{s=1}^t, h_t = i, (h_s)_{s=1}^t\\right),$$\n",
"каждая из которых интерпретируется как максимальная совместная вероятность появления известной наблюдаемой подпоследовательности $(x_s)_{s=1}^t$, появления произвольной скрытой подпоследовательности $(h_s)_{s=1}^t$ и события, заключающегося в том, что в $t$-й позиции скрытое состояние равнялось $i$.\n",
"\n",
"Также введём индексируемые позицией в последовательности $t$ и скрытым состоянием $i$ переменные\n",
"$$\\phi_{ti} = \\left(\\arg\\max_{(h_s)_{s=1}^t \\in \\{1, ..., m\\}^t}p\\!\\left((x_s)_{s=1}^t, h_t = i, (h_s)_{s=1}^t\\right)\\right)_{t-1},$$\n",
"каждая из которых интерпретируется как скрытое состояние в позиции $t-1$ в подпоследовательности скрытых состояний, максимизирующей выражение из правой части определения $\\psi_{ti}$.\n",
"\n",
"Алгоритм Витерби устроен так:\n",
"1. Для всех $i \\in \\{1, ..., m\\}$ вычислить начальные значения:\n",
"$$\\psi_{1i} = \\pi(i)B_{ix_1},$$\n",
"$$\\phi_{1i} = 0,$$\n",
"где 0 взят как произвольный некорректный номер состояния (номера состояний начинаются с 1). Иными словами, ни одну $\\phi_{1i}$ инициализировать не обязательно.\n",
"2. Для всех $(t, i) \\in \\{2, ..., T\\} \\times \\{1, ..., m\\} $ в порядке неубывания $t$ провести индуктивный шаг вперёд:\n",
"$$\\psi_{ti} = B_{ix_t} \\max_{j \\in \\{1, ..., m\\}}A_{ji}\\psi_{(t-1)j},$$\n",
"$$\\phi_{ti} = \\arg\\max_{j \\in \\{1, ..., m\\}}A_{ji}\\psi_{(t-1)j}.$$\n",
"3. Восстановить последнее скрытое состояние:\n",
"$$h_T = \\arg\\max_k \\psi_{Tk}.$$\n",
"4. Восстановить все остальные скрытые состояния индуктивными шагами назад:\n",
"$$h_t = \\phi_{(t+1)h_{t+1}}.$$\n",
"\n",
"Полученная алгоритмом Витерби последовательность $(h_t)_{t=1}^T$ является последовательностью с наибольшей совместной вероятностью $p\\!\\left((h_t)_{t=1}^T, (x_t)_{t=1}^T\\right)$ и, как следствие, с наибольшей условной вероятностью $p\\!\\left(\\left.(h_t)_{t=1}^T \\right| (x_t)_{t=1}^T\\right)$.\n",
"\n",
"Алгоритм Витерби делается устойчивым к проблеме округления малых чисел до нуля путём перехода к логарифмам, из-за которого произведения заменяются на суммы. Поскольку в выражениях, фигурирующих в описании алгоритма, операция сложения не участвует, а есть только операция умножения, более сложные приёмы (такие как, например, в алгоритме прохода вперёд) не нужны.\n",
"\n",
"#### Задача обучения скрытой марковской модели\n",
"\n",
"Есть несколько подходов к обучению скрытых марковских моделей:\n",
"* Алгоритм Баума-Уэлша, являющийся частным случаем [EM-алгоритма](__home_url__/notes/Общий взгляд на EM-алгоритм), применяемого во многих задачах обучения вероятностных графических моделей;\n",
"* Градиентный спуск в пространстве троек $(\\pi, A, B)$, после каждого шага которого строки матриц нормируются, чтобы их сумма, как это и требуется, равнялась 1.\n",
"\n",
"Помимо батчевых версий указанных алгоритмов обучения, сразу обучающихся на всём множестве имеющихся последовательностей без возможности дальнейшего дообучения, существуют и онлайновые версии, способные дообучаться на одной или более последовательностях.\n",
"\n",
"Что касается того, какой алгоритм обучения выбрать, оба алгоритма являются численными и оба не имеют гарантий того, что при оптимизации не произойдёт застревание около точки локального экстремума. По умолчанию алгоритм Баума-Уэлша предпочтительнее, однако встраивание в него регуляризации, необходимой для борьбы с присущим этому алгоритму переобучением, требует дополнительных математических изысканий, тогда как в случае с градиентным спуском достаточно просто добавить штрафное слагаемое к минимизируемой функции потерь.\n",
"\n",
"Для алгоритма Баума-Уэлша требуется в том числе делать проход назад при условии как-то зафиксированных параметров $\\pi$, $A$ и $B$. По аналогии с алгоритмом прохода вперёд введём следующие индексируемые позицией в последовательности $t$ и скрытым состоянием $i$ переменные:\n",
"$$\\beta_{ti} = p\\!\\left(\\left. (x_s)_{s=t+1}^T \\right| h_t = i \\right).$$\n",
"Вычисляются эти переменные индуктивными шагами назад от известного конца:\n",
"$$\\beta_{Ti} = 1,$$\n",
"$$\\beta_{ti} = \\sum_{j=1}^m A_{ij}B_{jx_{t+1}}\\beta_{(t+1)j}.$$\n",
"\n",
"Теперь можно описать, как устроен алгоритм Баума-Уэлша для обучения на множестве, состоящем ровно из одной последовательности.\n",
"\n",
"Сначала оценки параметров $\\pi$, $A$ и $B$ инициализируются случайно или на базе априорных знаний, а затем поочерёдно до сходимости проводятся следующие итерации EM-алгоритма:\n",
"* E-шаг:\n",
" - при текущих оценках $\\pi$, $A$ и $B$ вычислить все $\\alpha_{ti}$ алгоритмом прохода вперёд и все $\\beta_{ti}$ алгоритмом прохода назад;\n",
" - вычислить все условные вероятности того, что в позиции $t$ скрытое состояние равнялось $i$ при условии пронаблюдавшейся последовательности $(x_t)_{t=1}^T$ и при текущих оценках $\\pi$, $A$ и $B$:\n",
"$$\\gamma_{ti} = p\\!\\left(h_t = i \\, \\left| \\, (x_t)_{t=1}^T\\right.\\right) = \\frac{p\\!\\left(h_t = i, (x_t)_{t=1}^T\\right)}{p\\!\\left((x_t)_{t=1}^T\\right)} = \\frac{\\alpha_{ti}\\beta_{ti}}{\\sum_{j=1}^m \\alpha_{tj}\\beta_{tj}}.$$\n",
" - вычислить все условные вероятности перехода в позиции с $t$-й на $(t+1)$-ю из состояния $i$ в состояние $j$ при условии пронаблюдавшейся последовательности $(x_t)_{t=1}^T$ и при текущих оценках $\\pi$, $A$ и $B$:\n",
"$$\\xi_{tij} = p\\!\\left(h_t = i, h_{t+1} = j \\, \\left| \\, (x_t)_{t=1}^T\\right.\\right) = \\frac{p\\!\\left(h_t = i, h_{t+1} = j, (x_t)_{t=1}^T\\right)}{p\\!\\left((x_t)_{t=1}^T\\right)} = \\frac{\\alpha_{ti}A_{ij}\\beta_{(t+1)j}B_{jx_{t+1}}}{\\sum_{k=1}^m\\sum_{l=1}^m \\alpha_{tk}A_{kl}\\beta_{(t+1)l}B_{lx_{t+1}}}.$$\n",
"* M-шаг:\n",
" - Обновить $\\pi$ так:\n",
"$$\\pi(i) := \\gamma_{1i},$$\n",
"потому что справа стоит рассчитанная вероятность начать из состояния $i$;\n",
" - Обновить матрицу $A$ по следующей формуле обновления её элементов:\n",
"$$A_{ij} := \\frac{\\sum_{t=1}^{T-1} \\xi_{tij}}{\\sum_{t=1}^{T-1} \\gamma_{ti}},$$\n",
"потому что справа стоит отношение рассчитанного числа переходов из $i$ в $j$ к рассчитанному числу всех переходов из скрытого состояния $i$ (именно потому что состояние $i$ должно быть покинуто, сумма берётся лишь до $T-1$);\n",
" - Обновить матрицу $B$ по следующей формуле обновления её элементов:\n",
"$$B_{io} := \\frac{\\sum_{t=1}^T[x_t = o]\\,\\gamma_{ti}}{\\sum_{t=1}^T\\gamma_{ti}},$$\n",
"где квадратные скобки обозначают нотацию Айверсона, а справа стоит отношение рассчитанного числа нахождений в скрытом состоянии $i$, таких что было выпущено наблюдаемое значение $o$, к рассчитанному числу всех нахождений в скрытом состоянии $i$.\n",
"\n",
"Если же есть $N$ последовательностей в обучающей выборке, то:\n",
"* на E-шаге для всех $k$ от 1 до $N$ для $k$-й последовательности считаются свои значения $\\alpha_{ti}^{(k)}$, $\\beta_{ti}^{(k)}$, $\\gamma_{ti}^{(k)}$ и $\\xi_{tij}^{(k)}$, а также алгоритмом прохода вперёд вычисляется вероятность этой $k$-й последовательности, далее обозначаемая как $P_k$;\n",
"* на M-шаге вклад каждой последовательности учитывается с весом, обратно пропорциональным вероятности этой последовательности; например, формула обновления матрицы переходов между скрытыми состояниями будет выглядеть так:\n",
"$$A_{ij} := \\frac{\\sum_{k=1}^N \\left(\\frac{1}{P_k} \\sum_{t=1}^{T_k-1} \\xi_{tij}^{(k)}\\right)}{\\sum_{k=1}^N \\left(\\frac{1}{P_k}\\sum_{t=1}^{T_k-1} \\gamma_{ti}^{(k)}\\right)}.$$\n",
"Более сильная подгонка под те последовательности, которые кажутся маловероятными относительно текущих оценок $\\pi$, $A$ и $B$, позволяет обучать параметры так, чтобы все последовательности обучающей выборки выглядели вероятными."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"графические_модели",
"временные_ряды"
]
},
"source": [
"## Линейные динамические системы\n",
"\n",
"#### Введение\n",
"\n",
"Линейная динамическая система (ЛДС) — это вероятностная графическая модель, задаваемая сущностями из следующего списка:\n",
"* Последовательность $(h_t)_{t=1}^{\\infty}$ скрытых переменных, являющихся вещественными векторами, $h_t \\in \\mathbb{R}^m$;\n",
"* Последовательность $(x_t)_{t=1}^{\\infty}$ наблюдаемых переменных, являющихся вещественными векторами, $x_t \\in \\mathbb{R}^n$;\n",
"* Вероятностное распределение $p(h_t \\, \\vert \\, h_{t-1})$, описывающее переходы между скрытыми состояниями, причём предполагается, что это распределение является гауссовским и имеет вид $p(h_t \\, \\vert \\, h_{t-1}) = \\mathcal{N}_{Ah_{t-1}, \\, \\Gamma}(h_t)$, где $A$ и $\\Gamma$ — матрицы размера $m \\times m$;\n",
"* Вероятностное распределене $p(x_t \\, \\vert \\, h_t)$, задающее модель шума, из-за которого наблюдаемые значения отличаются от скрытых значений, причём предполагается, что это распределение тоже является гауссовским и имеет вид $p(x_t \\, \\vert \\, h_t) = \\mathcal{N}_{Ch_t, \\, \\Sigma}(x_t)$, где $C$ — матрица размера $n \\times m$, а $\\Sigma$ — матрица размера $n \\times n$;\n",
"* Начальное распределение скрытого состояния $p(h_1) = \\mathcal{N}_{\\lambda_0, \\, W_0}(h_1)$, где $\\lambda_0$ — вектор-столбец размера $m \\times 1$, а $W_0$ — начальная ковариационная матрица размера $m \\times m$.\n",
"\n",
"С учётом всего вышесказанного вероятностная модель ЛДС принимает вид:\n",
"$$p\\left((h_t)_{t=1}^T, (x_t)_{t=1}^T \\, \\vert \\, A, \\Gamma, C, \\Sigma, \\lambda_0, W_0 \\right) = \\mathcal{N}_{\\lambda_0, \\, W_0}(h_1) \\mathcal{N}_{Ch_1, \\, \\Sigma}(x_1) \\prod_{t=2}^T \\left(\\mathcal{N}_{Ah_{t-1}, \\, \\Gamma}(h_t) \\mathcal{N}_{Ch_t, \\, \\Sigma}(x_t)\\right),$$\n",
"где $T$ — некоторый момент времени. Если расписать все плотности гауссовских распределений, то выкладкой можно получить, что совместное распределение последовательностей $(h_t)_{t=1}^T$ и $(x_t)_{t=1}^T$ и само является гауссовским. Как следствие, его маржинальные и условные распределения тоже являются гауссовскими.\n",
"\n",
"По сути, линейные динамические системы можно считать аналогами [скрытых марковских моделей](__home_url__/notes/Скрытые марковские модели) (СММ) для того случая, когда скрытые переменные являются непрерывными, а не дискретными. Что для ЛДС, что для СММ граф взаимозависимости переменных выглядит одинаково.\n",
"\n",
"Для ЛДС можно сформулировать четыре следующих задачи:\n",
"* Зная параметры $A$, $\\Gamma$, $C$, $\\Sigma$, $\\lambda_0$ и $W_0$, провести восстановление скрытого состояния в режиме реального времени, то есть научиться вычислять условные вероятности $p(h_t \\, \\vert \\, (x_i)_{i=1}^{t-1})$ и $p(h_t \\, \\vert \\, (x_i)_{i=1}^{t})$, где первую условную вероятность можно использовать для прогнозирования $t$-го скрытого состояния на $(t-1)$-м шаге, а вторую — для уточнения такого прогноза, после того как поступило $t$-e наблюдаемое значение.\n",
"* Зная параметры $A$, $\\Gamma$, $C$, $\\Sigma$, $\\lambda_0$ и $W_0$, найти наиболее вероятную последовательность скрытых состояний $(h_t)_{t=1}^T$ по заранее известной последовательности наблюдаемых состояний $(x_t)_{t=1}^T$.\n",
"* Имея последовательность наблюдаемых состояний $(x_t)_{t=1}^T$ и последовательность скрытых состояний $(h_t)_{t=1}^T$ обучить с учителем параметры $A$, $\\Gamma$, $C$, $\\Sigma$, $\\lambda_0$ и $W_0$, то есть найти их значения, максимизирующие полное правдоподобие.\n",
"* Имея последовательность наблюдаемых состояний $(x_t)_{t=1}^T$ обучить без учителя параметры $A$, $\\Gamma$, $C$, $\\Sigma$, $\\lambda_0$ и $W_0$, то есть найти их значения, максимизирующие неполное правдоподобие.\n",
"\n",
"#### Задача восстановления скрытого состояния в режиме реального времени\n",
"\n",
"Эту задачу будем решать индуктивно.\n",
"\n",
"На нулевом шаге прогноз для шага $t=1$ дан в виде начального распределения скрытого состояния: $p(h_1) = \\mathcal{N}_{\\lambda_0, \\, W_0}(h_1)$. Коррекция этого прогноза после поступления значения $x_1$, а именно $p(h_1 \\, \\vert \\, x_1)$, является некоторым гауссовским распределением в силу того, что, по определению, это условное матожидание гауссовской величины при условии другой гауссовской величины, получаемой из первой добавлением гауссовского шума. Матожидание и ковариационную матрицу для $p(h_1 \\, \\vert \\, x_1)$ можно найти по известным формулам для такого случая:\n",
"$$p(h_1 \\, \\vert \\, x_1) = \\mathcal{N}_{\\mu_1, \\, V_1}(h_1),$$\n",
"$$V_1 = (W_0^{-1} + C^T\\Sigma^{-1}C)^{-1},$$\n",
"$$\\mu_1 = V_1(C^T\\Sigma^{-1}x_1 + W_0^{-1}\\lambda_0),$$\n",
"где символ $T$ в верхнем регистре обозначает транспонирование, а не возведение в степень, равную некоторому моменту времени.\n",
"\n",
"Окончив с базой индукции, перейдём к шагу индукции. Предположим, что известно распределение $p(h_{t-1} \\, \\vert \\, (x_i)_{i=1}^{t-1})$, и по нему выведем распределение $p(h_t \\, \\vert \\, (x_i)_{i=1}^{t-1})$. Для начала отметим, что все такие распределения являются гауссовскими, то есть допускают представление в виде $p(h_{t} \\, \\vert \\, (x_i)_{i=1}^{t}) = \\mathcal{N}_{\\mu_{t}, \\, V_{t}}(h_{t})$, где через $\\mu_{t}$ и $V_{t}$ обозначены матожидание и ковариационная матрица.\n",
"\n",
"Вычисление прогноза скрытого состояния на следующий шаг выглядит так:\n",
"$$p(h_t \\, \\vert \\, (x_i)_{i=1}^{t-1}) = \\int p(h_t, h_{t-1} \\, \\vert \\, (x_i)_{i=1}^{t-1}) \\, dh_{t-1} = \\int p(h_t \\, \\vert \\, h_{t-1}) p(h_{t-1} \\, \\vert \\, (x_i)_{i=1}^{t-1}) \\, dh_{t-1} = \\int \\mathcal{N}_{Ah_{t-1}, \\, \\Gamma}(h_t) \\mathcal{N}_{\\mu_{t-1}, \\, V_{t-1}}(h_{t-1}) \\, dh_{t-1} = \\mathcal{N}_{A\\mu_{t-1}, \\, \\Gamma + AV_{t-1}A^T}(h_{t}).$$\n",
"\n",
"Глядя на получившуюся правую часть, введём пару новых обозначений:\n",
"$$\\overline{\\mu}_t = A\\mu_{t-1},$$\n",
"$$\\overline{V}_t = \\Gamma + AV_{t-1}A^T.$$\n",
"\n",
"Тем самым через $\\mu_{t-1}$, $V_{t-1}$ и параметры ЛДС, предполагающиеся известными, удалось выразить $\\overline{\\mu}_t$ и $\\overline{V}_t$. Они отличаются от нужных для проведения шага индукции $\\mu_t$ и $V_t$ тем, что в них не учитывается информация о текущем наблюдаемом значении $x_t$ (в частности, можно считать, что $\\lambda_0 = \\overline{\\mu}_1$ и $W_0 = \\overline{V}_1$). Чтобы выразить $\\mu_t$ и $V_t$ через $\\overline{\\mu}_t$ и $\\overline{V}_t$, проведём коррекцию, вносимую получением информации об $x_t$.\n",
"\n",
"Для данной коррекции введём ещё одно обозначение, сокращающее размер формул:\n",
"$$K_t = \\overline{V}_tC^T(C\\overline{V}_tC^T + \\Sigma)^{-1}.$$\n",
"\n",
"С использованием этого обозначения можно получить такие формулы коррекции:\n",
"$$\\mu_t = \\overline{\\mu}_t + K_t(x_t - C\\overline{\\mu}_t),$$\n",
"$$V_t = (I - K_tC)\\overline{V}_t,$$\n",
"где $I$ — единичная матрица размера $m \\times m$. Вывод этих формул опущен из-за громоздкости.\n",
"\n",
"Напомним, что начальные значения $\\mu_1$ и $V_1$ были вычислены при описании базы индукции, так что теперь можно вычислить $\\mu_t$ и $V_t$ для любого $t$. Для этого, начав с прогноза на шаг $t=1$, можно в цикле сначала прогнозировать на следующий шаг, а потом корректировать такой прогноз на соответствующее наблюдаемое значение.\n",
"\n",
"Подобная процедура называется проходом вперёд (forward algorithm) — так же, как алгоритм для скрытых марковских моделей, который вычислял аналогичные распределения (правда, с другой конечной целью: с целью вычислительно эффективно найти вероятность последовательности наблюдаемых состояний). Также данный алгоритм называют фильтром Калмана. Слово «фильтр» здесь интерпретируется в том ключе, что наблюдаемые значения предполагаются зашумлёнными версиями истинных скрытых значений и требуется сгладить временной ряд, отфильтровав шум.\n",
"\n",
"#### Задача нахождения наиболее вероятной последовательности скрытых состояний\n",
"\n",
"В отличие от скрытых марковских моделей для решения этой задачи отдельный алгоритм (в случае СММ называющийся алгоритмом Витерби) не нужен: достаточно алгоритма прохода вперёд-назад. Для ЛДС наиболее вероятная последовательность скрытых значений является последовательностью, составленной из индивидуально наиболее вероятных скрытых значений. Это так, потому что, как отмечалось выше, совместное распределение $p\\left((h_t)_{t=1}^T, (x_t)_{t=1}^T \\, \\vert \\, A, \\Gamma, C, \\Sigma, \\lambda_0, W_0 \\right)$ является нормальным, а для многомерного нормального распределения проекции моды на ось какой-либо компоненты совпадает с модой проекции распределения на эту ось (можно ещё добавить, что мода совпадает с матожиданием).\n",
"\n",
"Тем самым задача сводится к такой: по известной последовательности $(x_i)_{i=1}^T$ найти моду распределения $p(h_t \\, \\vert \\, (x_i)_{i=1}^T)$, где $1 \\le t \\le T$. Иными словами, теперь требуется восстановить скрытое состояние, зная не только прошлые и настоящие наблюдения, но и некоторые последующие наблюдения. Интуитивно говоря, проход назад для того и нужен, чтобы учесть информацию о более поздних наблюдаемых значениях.\n",
"\n",
"Алгоритм прохода вперёд-назад состоит из двух этапов:\n",
"* Сначала алгоритмом прохода вперёд, разобранным ранее, вычисляются вектор среднего и ковариационная матрица для распределения $p(h_T \\, \\vert \\, (x_i)_{i=1}^{T}) = \\mathcal{N}_{\\mu_T, \\, W_T}(h_T)$.\n",
"* Затем алгоритмом прохода назад, который будет разобран ниже, вычисляются векторы среднего и ковариационные матрицы для распределений $p(h_t \\, \\vert \\, (x_i)_{i=1}^{T})$ в порядке от $t = T-1$ к $t=1$.\n",
"\n",
"Введём обозначения $\\tilde{\\mu}_t$ и $\\tilde{V}_t$: $p(h_t \\, \\vert \\, (x_i)_{i=1}^{T}) = \\mathcal{N}_{\\tilde{\\mu}_t, \\, \\tilde{V}_t}(h_t)$. В соответствии с этими обозначениями $\\tilde{\\mu}_T = \\mu_T$ и $\\tilde{V}_T = V_T$, а для остальных $t$ эти параметры ищутся при проходе назад.\n",
"\n",
"Для $t$, лежащего между $T-1$ и 1 шаг в алгоритме прохода назад выглядит так:\n",
"$$J_t = V_t A^T \\overline{V}_t^{-1},$$\n",
"$$\\tilde{\\mu}_t = \\mu_t + J_t (\\tilde{\\mu}_{t+1} - A\\mu_t),$$\n",
"$$\\tilde{V}_t = V_t + J_t (\\tilde{V}_{t+1} - \\overline{V}_t)J_t^T,$$\n",
"где $J_t$ — вспомогательное обозначение, сокращающее формулы, а $V_t$ и $\\overline{V}_t$ были вычислены при проходе вперёд. Выписанную тройку формул называют РТС-уравнениями по первым буквам имён её авторов."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"генеративные_модели",
"EM"
]
},
"source": [
"## Восстановление плотности при помощи EM-алгоритма\n",
"\n",
"Пусть дана обучающая выборка из $l$ независимых одинаково распределённых объектов, представленная матрицей $X$ размера $l \\times n$, где $n$ — количество признаков. Пусть также имеется предположение об аналитическом виде функции плотности распределения, откуда порождена выборка, то есть имеется $p(x \\, \\vert \\, \\theta)$, где $x$ — вектор-строка признакового описания объекта, а $\\theta$ — вектор неизвестных параметров распределения. Например, вид такой функции может быть известен благодаря представлениям о предметной области.\n",
"\n",
"Если найти $\\theta$, то тем самым будет восстановлена плотность, что, с одной стороны, позволит углубить понимание данных и, с другой стороны, даст возможность сэмплировать новые объекты из того же самого распределения.\n",
"\n",
"Одним из возможных решений описанной задачи восстановления плотности является метод максимального правдоподобия. Впрочем, если распределение $p(x \\, \\vert \\, \\theta)$ не принадлежит к классу экспоненциальных распределений, могут возникнуть трудности, связанные с многоэкстремальностью возникающей оптимизационной задачи.\n",
"\n",
"Есть и более интересное решение. Можно дополнительно ввести $h$ скрытых признаков объекта. Обозначим вектор-строку скрытых признаков объекта за $t$ и переопределим $\\theta$ как вектор параметров следующей модели:\n",
"$$p(x, t \\, \\vert \\, \\theta) = p(x \\, \\vert \\, t, \\theta) p(t \\, \\vert \\, \\theta).$$\n",
"Если сделать это удачным образом, то можно добиться того, что распределение $p(x, t \\, \\vert \\, \\theta)$ окажется в экспоненциальном классе. Правда, из-за присутствия скрытых переменных найти $\\theta$ методом максимального правдоподобия уже не получится. Придётся использовать [EM-алгоритм](__home_url__/notes/Общий взгляд на EM-алгоритм).\n",
"\n",
"Выше было употреблено размытое словосочетание «удачным образом». Проиллюстрируем его значение на примерах. Допустим, исходная $p(x \\, \\vert \\, \\theta)$ соответствует смеси нормальных распределений. Это распределение не принадлежит к экспоненциальному классу. Однако предположим, что теперь введён один дополнительный категориальный признак $t$, кодирующий номер компоненты смеси, откуда был порождён объект. Тогда $p(x, t \\, \\vert \\, \\theta)$ будет лежать в экспоненциальном классе, потому что разлагается в произведение двух множителей, из которых один, а именно $p(x \\, \\vert \\, t, \\theta)$, является плотностью нормального распределения, а другой, то есть $p(t \\, \\vert \\, \\theta)$, является плотностью дискретного распределения. Сам по себе такой шаг хорош. Однако и его можно реализовать неправильно. Количество компонент смеси $k$ должно быть подобрано, и здесь есть опасность недоподгонки или переподгонки. Например, в вырожденном случае $k = l$ каждый объект окажется в своей собственной компоненте, ковариационную матрицу которой, кстати, нельзя будет достоверно оценить."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"гауссовские_процессы",
"байесовские_методы",
"непараметрические_методы",
"пререквизиты_из_математики"
]
},
"source": [
"## Регрессия на основе гауссовских процессов\n",
"\n",
"#### Необходимые определения\n",
"\n",
"Регрессия на основе гауссовских процессов опирается на предположение о том, что целевая переменная порождена гауссовским процессом. Значит, уместно дать определение того, что такое гауссовский процесс. Сделаем это поэтапно.\n",
"\n",
"Определим случайное поле как множество мощности континуум, состоящее из случайных величин $\\xi(x, \\omega)$, где $x \\in \\mathbb{R}^n$ индексирует множество, $n \\in \\mathbb{N}$, а $\\omega$ обозначает тот или иной исход случайности и далее ради краткости будет опускаться. В контексте машинного обучения каждый индексирующий элемент $x$ интерпретируется как признаковое описание некого объекта, и тогда $n$ становится количеством признаков. Если считать, что целевая переменная задаётся случайным полем, то для некоторого объекта $x$ его пронаблюдавшаяся целевая переменная является реализацией случайной величины $\\xi(x, \\omega)$ при выпавшем исходе $\\omega$. Если же зафиксировать $\\omega$ на исходе $\\overline{\\omega}$, то тогда $\\xi(x, \\overline{\\omega})$ будет детерминированной функцией, описывающей зависимость реализовавшейся целевой переменной от признаков объекта. \n",
"\n",
"Гауссовским процессом назовём случайное поле, такое что $\\forall k \\in \\mathbb{N}$ $\\forall x_1, \\dots, x_k$, таких что все эти $x_i \\in \\mathbb{R}^n$, $k$-мерная случайная величина $(\\xi(x_1), \\dots, \\xi(x_k))$ имеет $k$-мерное гауссовское распределение. На самом деле, иногда гауссовские процессы называют гауссовскими полями, и такой термин звучит более корректно, потому что отсылает к пространственной протяжённости (к изменению каждого из признаков), а не к протяжённости во времени.\n",
"\n",
"Подобно тому, как гауссовское распределение полностью задаётся математическим ожиданием и ковариационной матрицей, гауссовский процесс полностью задаётся парой из двух функций, из которых первая называется функцией среднего, а вторая — ковариационной функцией:\n",
"$$m(x) = \\mathbb{E}\\xi(x),$$\n",
"$$K(x_1, x_2) = \\mathrm{Cov}(\\xi(x_1), \\xi(x_2)).$$\n",
"\n",
"Важным классом гауссовских процессов являются стационарные процессы. Стационарностью (в слабом смысле) называют соблюдение следующих двух свойств: $m(x)$ не зависит от $x$, а $K(x_1, x_2)$ зависит только от $x_1 - x_2$, то есть представима в виде $K(x_1 - x_2)$.\n",
"\n",
"#### Формулы Андерсона\n",
"\n",
"Пусть целевая переменная порождена стационарным гауссовским процессом $\\xi(x)$ с $m(x) \\equiv 0$. Предположение о стационарности введено для упрощения выкладок, и оно не ограничивает общность, потому что произвольный гауссовский процесс можно остационарить. Опять же, из стационарного процесса всегда можно вычесть его среднее, чтобы получить процесс с нулевым средним.\n",
"\n",
"Далее, предположим, что ковариционная функция $K(x_1 - x_2)$ известна в аналитическом виде. На практике эту функцию можно просто выбрать, и такой выбор будет в чём-то похож на выбор ширины окна при использовании метода окна Парзена.\n",
"\n",
"Задача заключается в том, что, имея обучающую выборку $x_1, \\dots, x_l$, для которой известны целевые переменные $\\xi(x_1), \\dots, \\xi(x_l)$, требуется для новой выборки $x_{l+1}, \\dots, x_s$ при сделанных предположениях найти условное распределение $p(\\xi(x_{l+1}), \\dots, \\xi(x_s) \\, \\vert \\, \\xi(x_1), \\dots, \\xi(x_l))$.\n",
"\n",
"Ничего сложного в такой задаче нет:\n",
"* искомая условная плотность записывается как отношение совместной плотности к плотности условия;\n",
"* эти две плотности, по определению гауссовского процесса, являются плотностями $s$-мерного и $l$-мерного гауссовских распределений;\n",
"* отношением двух гауссиан является гауссиана, параметры который находятся арифметической выкладкой.\n",
"\n",
"Если проделать все выкладки, получится, что распределение $p(\\xi(x_{l+1}), \\dots, \\xi(x_s) \\, \\vert \\, \\xi(x_1), \\dots, \\xi(x_l))$ является $(s - l)$-мерным гауссовским распределением с вектором среднего $S^TC^{-1}\\overrightarrow{\\xi}$ и ковариационной матрицей $Q - S^TC^{-1}S$, где:\n",
"* $C$ — матрица размера $l \\times l$, элементы которой определены как $c_{ij} = K(x_i - x_j)$, то есть $C$ является ковариационной матрицей проекции гауссовского процесса на объекты обучающей выборки;\n",
"* $S$ — матрица размера $l\\times(s-l)$, элементы которой определены как $s_{ij} = K(x_i - x_{l+j})$, то есть $S$ отражает ковариацию между объектами обучающей выборки и новыми объектами;\n",
"* $Q$ — матрица размера $(s-l)\\times(s-l)$, которой определены как $q_{ij} = K(x_{l+i} - x_{l+j})$, то есть $Q$ является ковариационной матрицей проекции гауссовского процесса на новые объекты;\n",
"* $\\overrightarrow{\\xi}$ — вектор-столбец значений целевой переменной на объектах обучающей выборки.\n",
"\n",
"Формулы для среднего и ковариационной матрицы условного гауссовского процесса при условии обучающей выборки называют формулами Андерсона.\n",
"\n",
"#### Частные случаи\n",
"\n",
"Если выбрана ковариационная функция $K(x_1 - x_2)$, являющаяся непрерывной в 0, то обусловленный гауссовский процесс, восстановленный по точкам обучающей выборки, будет иметь непрерывную функцию среднего, проходящую через точки обучающей выборки, а дисперсии случайных величин, соответствующих точкам обучающей выборки, будут нулевыми: $\\mathrm{Var} \\, \\xi(x_i) = 0$ $\\forall i \\in \\{1, \\dots, l\\}$. Это следует из формул Андерсона.\n",
"\n",
"Если же $K(x_1 - x_2)$ непрерывна всюду кроме 0, где имеет устранимый разрыв (то есть имеет предел слева при стремлении к 0, имеет предел справа при стремлении к 0, и эти два предела равны), то обусловленный гауссовский процесс будет иметь разрывную функцию среднего, имеющую устранимые разрывы в точках обучающей выборки и проходящую через точки обучающей выборки благодаря этим разрывам. Любая функция $K(x_1 - x_2)$ с приведёнными выше свойствми является суммой некоторой непрерывной ковариационной функции и ковариационной функции белого шума, так что исходный гауссовский процесс можно рассматривать как гауссовский процесс с зашумлением. Таким образом, для обусловленной $m(x)$ пределы $\\lim_{x \\rightarrow x_i} m(x)$ могут пониматься как очищенные от шума значения целевой переменной в точках обучающей выборки.\n",
"\n",
"#### Выбор ковариационной функции через принцип наибольшей обоснованности\n",
"\n",
"Допустим, ковариационная функция ищется в каком-либо параметрическом семействе. Часто используют следующее семейство: $$K(x_1 - x_2) = Ae^{-\\gamma\\Vert x_1 - x_2\\Vert^2} + \\sigma[x_1 = x_2],$$\n",
"где запись $\\Vert\\cdot\\Vert$ обозначает евклидово расстояние, а квадратные скобки обозначают нотацию Айверсона. Параметры этого семейства интерпретируются так:\n",
"* $A$ задаёт общий уровень ковариации;\n",
"* $\\gamma$ задаёт степень затухания ковариации;\n",
"* $\\sigma$ задаёт степень зашумлённости, ведь слагаемое с этим параметром является ковариационной функцией белого шума.\n",
"\n",
"Возникает вопрос, как подобрать параметры ковариационной функции. Ответ на него может быть дан байесовским [принципом наибольшей обоснованности](__home_url__/notes/Принцип наибольшей обоснованности). Согласно этому принципу, следует выбрать ту тройку $(A, \\gamma, \\sigma)$, для которой максимальна обоснованность, то есть вероятность обучающей выборки. Поскольку исходный гауссовский процесс является стационарным с нулевым средним, обоснованность равна плотности в точке $(\\xi(x_1), \\dots, \\xi(x_l))$ $l$-мерного нормального распределения с нулевым вектором среднего и ковариационной матрицей $C$. Решение соответствующей оптимизационной задачи даёт рекомендованные параметры ковариационной функции.\n",
"\n",
"#### Непараметричность\n",
"\n",
"Даже если ковариационная функция зафиксирована заранее или ищется в параметрическом семействе, регрессия на основе гауссовских процессов является непараметрическим методом.\n",
"\n",
"Здесь предлагается посмотреть на непараметричность не в терминах того, что множество зависимостей, могущих быть выраженными алгоритмом обучения, всюду плотно во множестве всех зависимостей, а в терминах того, что сложность восстановленной зависимости определяется не только гиперпараметрами, но и размером выборки.\n",
"\n",
"Возьмём, скажем, случайный лес, у которого зафиксировали количество деревьев на уровне $q$ и ограничили максимальную глубину каждого из деревьев уровнем $d$. Такой случайный лес заведомо не сможет выдавать более $2^{qd}$ различных предсказаний. Однако если обучать регрессию на основе гауссовских процессов, чем больше обучающая выборка, тем более сложным может быть обусловленный процесс. Например, если взять какую-то одну непрерывную ковариационную функцию и зафиксировать её, предсказания будут содержать не меньше уникальных значений, чем их было у целевой переменной в обучающей выборке (по той хотя бы причине, что на объектах обучающей выборки предсказывается сама их целевая переменная)."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"пререквизиты_из_математики"
]
},
"source": [
"## Распределение Дирихле и случайные процессы Дирихле\n",
"\n",
"#### Определение и свойства распределения Дирихле\n",
"\n",
"Для $k \\in \\mathbb{N}$ рассмотрим множество $P$ векторов $\\pi$ длины $k$, таких что $\\sum_{i=1}^k \\pi_i = 1$ и $\\forall i: 1 \\le i \\le k,$ соблюдается неравенство $\\pi_i \\ge 0$. Любой такой вектор $\\pi \\in P$ может интерпретироваться как вектор вероятностей, которым параметризуется категориальное распределение с $k$ исходами.\n",
"\n",
"На множестве $P$ и определено распределение Дирихле. Это распределение параметризуется вектором $\\alpha$ длины $k$, таким что $\\forall i: 1 \\le i \\le k,$ верно неравенство $\\alpha_i > 0$. Обозначается распределение Дирихле как $\\mathrm{Dir}_{\\alpha}(\\pi)$, и его плотность в точке $\\pi$ имеет следующий вид:\n",
"$$\\mathrm{Dir}_{\\alpha}(\\pi) = \\frac{\\Gamma(\\sum_{i=1}^k \\alpha_i)}{\\prod_{i=1}^k \\Gamma(\\alpha_i)}\\prod_{i=1}^k\\pi_i^{\\alpha_i - 1},$$\n",
"где $\\Gamma$ — гамма-функция.\n",
"\n",
"Первый множитель правой части является нормировочной константой, и поведение функции плотности распределения Дирихле определяется множителем $\\prod_{i=1}^k\\pi_i^{\\alpha_i - 1}$. Из вида этого множителя вытекает, что если все $\\alpha_i = 1$, то распределение Дирихле является равномерным распределением на множестве $P$, если все $\\alpha_i > 1$, то распределение Дирихле имеет плотность, сконцентрированную около векторов $\\pi$, таких что одна их компонента равна 1, а остальные — 0. Наконец, если все $\\alpha_i < 1$, плотность распределения Дирихле имеет моду в векторе $\\pi$, таком что в соответствующем ему категориальном распределении все исходы более-менее вероятны.\n",
"\n",
"Из свойств распределения Дирихле можно отметить следующие:\n",
"* Агрегационное свойство: если $\\pi \\sim \\mathrm{Dir}_{\\alpha}(\\pi)$, то $(k-1)$-мерный вектор $(\\pi_1 + \\pi_2, \\pi_3, \\dots, \\pi_k)$ имеет распределение Дирихле, параметризованное $(k-1)$-мерным вектором $(\\alpha_1 + \\alpha_2, \\alpha_3, \\dots, \\alpha_k)$:\n",
"$$(\\pi_1 + \\pi_2, \\pi_3, \\dots, \\pi_k) \\sim \\mathrm{Dir}_{(\\alpha_1 + \\alpha_2, \\alpha_3, \\dots, \\alpha_k)}((\\pi_1 + \\pi_2, \\pi_3, \\dots, \\pi_k)).$$\n",
"Предельный случай этого свойства принимает вид:\n",
"$$(\\pi_1, \\sum_{i=2}^k \\pi_i) \\sim \\mathrm{Beta}_{\\alpha_1, \\sum_{i=2}^k \\alpha_i}(\\pi_1) = \\frac{1}{B(\\alpha_1, \\sum_{i=2}^k \\alpha_i)}\\pi_1^{\\alpha_1 - 1}(1 - \\pi_1)^{\\sum_{i=2}^k \\alpha_i - 1},$$\n",
"где $\\mathrm{Beta}$ — бета-распределение (его плотность выписана после знака равенства), а $B$ — бета-функция. Отметим, что бета-распределение является двумерным аналогом распределения Дирихле.\n",
"* Свойство условной независимости: если $\\pi \\sim \\mathrm{Dir}_{\\alpha}(\\pi)$, то его условное распределение при условии одной компоненты по-прежнему является распределением Дирихле:\n",
"$$p(\\pi_2, \\dots, \\pi_k \\, \\vert \\, \\pi_1) = \\frac{p(\\pi)}{p(\\pi_1)} = \\frac{\\mathrm{Dir}_{\\alpha}(\\pi)}{ \\mathrm{Beta}_{\\alpha_1, \\sum_{i=2}^k \\alpha_i}(\\pi_1)} = \\frac{1}{\\mathrm{const}}\\prod_{i=1}^k\\left(\\frac{\\pi_i}{1 - \\pi_1}\\right)^{\\alpha_i - 1} = \\mathrm{Dir}_{(\\alpha_2, \\dots, \\alpha_k)}\\left((\\frac{\\pi_2}{1 - \\pi_1}, \\dots, \\frac{\\pi_k}{1 - \\pi_1})\\right).$$\n",
"\n",
"#### Сэмплирование из распределения Дирихле\n",
"\n",
"Чтобы сэмплировать векторы $\\pi$ из распределения Дирихле с параметром $\\alpha$, можно сгенерировать $k$ величин из гамма-распределения, а потом их отнормировать. Однако в контексте случайных процессов Дирихле более интересен другой способ, называющийся stick breaking.\n",
"\n",
"Устроен этот способ так. Сначала сэмплируются $k$ величин $v_1, \\dots, v_k$ из следующих бета-распределений:\n",
"$$v_1 \\sim \\mathrm{Beta}_{\\alpha_1, \\sum_{i=2}^k\\alpha_i}(v_1),$$\n",
"$$v_2 \\sim \\mathrm{Beta}_{\\alpha_2, \\sum_{i=3}^k\\alpha_i}(v_2),$$\n",
"$$\\dots,$$\n",
"$$v_{k-1} \\sim \\mathrm{Beta}_{\\alpha_{k-1}, \\alpha_k}(v_{k-1}),$$\n",
"$$v_{k} = 1.$$\n",
"Далее $i$-я компонента сэмплируемого вектора $\\pi$ находится по формуле:\n",
"$$\\pi_i = \\prod_{j=1}^{i-1}(1-v_j) v_i.$$\n",
"\n",
"Корректность работы такого способа доказывается через агрегационное свойство распределения Дирихле. Интуитивно говоря, с каждым $v_i$ определяется, какая доля ещё нераспределённой вероятности достанется $i$-му исходу.\n",
"\n",
"#### Определение и свойства случайного процесса Дирихле\n",
"\n",
"Ниже даётся неконструктивное определение процесса Дирихле. Иными словами, это определение не позволяет построить сам процесс, из-за чего для такого определения доказана сложная теорема, гласящая, что случайный процесс с перечисленными свойствами существует и единственен.\n",
"\n",
"Итак, случайный процесс Дирихле — это случайный процесс $\\xi$, задаваемый следующими свойствами:\n",
"* индексирующим элементом континуального множества случайных величин является множество $A$, являющееся подмножеством некоторого множества $U$. Для любого $A \\subset U$ определена случайная величина $\\xi(A, \\omega)$, задающая значение процесса на $A$ при исходе $\\omega$ (далее зависимость от исхода иногда будет опускаться). Поскольку процесс индексируется подмножеством, а не одномерным действительным числом, интерпретируемым как время, иногда его называют случайным полем Дирихле;\n",
"* каким бы ни оказался зафиксированный исход $\\overline{\\omega}$, реализация процесса Дирихле при этом исходе будет вероятностной мерой на $U$. Иными словами, функция $A \\mapsto \\xi(A, \\overline{\\omega})$ является вероятностной мерой на $U$;\n",
"* случайный процесс Дирихле параметризуется следующими двумя параметрами:\n",
" - $G$ — вероятностная мера на $U$, интерпретируемая как аналог среднего,\n",
" - $\\alpha > 0$ — неотрицательное число, контролирующее разреженность;\n",
"* одномерной проекцией процесса Дирихле на некоторое зафиксированное подмножество $\\overline{A}$ является следующая случайная величина:\n",
"$$\\xi(\\overline{A}) \\sim \\mathrm{Beta}_{\\alpha G(\\overline{A}), \\alpha(1 - G(\\overline{A}))}(\\xi(\\overline{A})),$$\n",
"откуда видно, что $G$ интерпретируется как среднее, потому что математическое ожидание $\\xi(\\overline{A})$ равно $G(\\overline{A})$;\n",
"* если взять $m$ подмножеств $A_1, \\dots, A_m$, таких что $\\bigcup_{i=1}^m A_i = U$ и $\\forall i, j: 1 \\le i < j \\le m$, $A_i \\cap A_j = \\emptyset$, то $m$-мерной проекцией процесса Дирихле на $(A_1, \\dots, A_m)$ будет $m$-мерная случайная величина из распределения Дирихле:\n",
"$$(\\xi(A_1), \\dots, \\xi(A_m)) \\sim \\mathrm{Dir}_{(\\alpha G(A_1), \\dots, \\alpha G(A_m))}((\\xi(A_1), \\dots, \\xi(A_m))).$$\n",
"\n",
"Из определения процесса Дирихле вытекает одно его важное свойство: мера $\\xi(A)$, являющаяся некоторой реализацией, не может иметь непрерывной части. Доказывается это от противного. Если бы у этой меры была непрерывная часть, можно бы было подобрать два подмножества $A_1$ и $A_2$ так, чтобы знание $\\xi(A_1)$ давало бы оценку на $\\xi(A_2)$, а это противоречило бы свойству условной независимости распределения Дирихле. Оставив в стороне существование сингулярных мер, будем считать, что отсюда выводится то, что мера $\\xi(A)$ является атомарной, то есть вся вероятность сосредоточена в некотором счётном множестве элементов $U$. Это свойство процесса Дирихле позволяет сэмплировать из него.\n",
"\n",
"#### Сэмплирование из случайного процесса Дирихле\n",
"\n",
"Сэмплирование из процесса Дирихле — получение различных реализаций процесса Дирихле при каких-либо исходах. Как было показано в предыдущем разделе, любую реализацию процесса Дирихле можно представить в виде счётного множества пар $\\{(\\theta_i, \\pi_i)\\}_{i=1}^{+\\infty}$, где $\\theta_i \\in U$, а $\\pi_i$ образуют вероятности, то есть $\\sum_{i=1}^{+\\infty}\\pi_i = 1$ и $\\pi_i > 0$. Поскольку на практике невозможно ждать бесконечно долго, критерием остановки порождения какой-либо реализации является получение множества пар такой длины, что сумма его $\\pi_i$ отличается от 1 менее чем на пределы машинной точности.\n",
"\n",
"Существует два наиболее известных способа сэмплировать из случайного процесса Дирихле:\n",
"* схема китайских ресторанов;\n",
"* stick breaking.\n",
"\n",
"Схема китайских ресторанов устроена следующим образом:\n",
"* $\\theta_1 \\in U$ сэплируется из $G$, этому элементу приписывается счётчик $n_1$, инициализированный 1, а вероятность этого элемента $\\pi_1$ будет вычисляться по формуле $\\pi_1 = \\frac{n_1}{n_1 + \\alpha}$.\n",
"* Целочисленная переменная $k_*$ инициализируется 2; интерпретируется эта переменная как количество имеющихся на текущем этапе различных $\\theta_i$, увеличенное на 1;\n",
"* Далее вплоть до достижения условия остановки вида $\\frac{\\alpha}{\\sum_{i=1}^{k_* - 1} n_i + \\alpha} < \\varepsilon$, где $\\varepsilon$ — порог точности, происходит следующее:\n",
" - с вероятностью $\\pi_i = \\frac{n_i}{n_i + \\alpha}$ может быть выбран один из ранее просэмплированных $\\theta_i$, тогда его счётчик $n_i$ увеличивается на 1 (что повышает вероятность выбора этого элемента в дальнейшем);\n",
" - с вероятностью $\\frac{\\alpha}{\\sum_{i=1}^{k_* - 1} n_i + \\alpha}$ может быть выбрано сэмплирование нового элемента, тогда $\\theta_{k_*}$ сэмплируется из $G$, ему приписывается счётчик $n_{k_*}$, инициализированный 1, а сама переменная $k_*$ увеличивается на 1, чтобы обозначать уже потенциальный следующий элемент.\n",
" \n",
"Схема stick breaking напоминает аналогичную схему для распределения Дирихле:\n",
"* счётчик $i$ инициализируется 1;\n",
"* пока $i = 1$ или не стало верным, что $\\prod_{j=1}^{i-1}(1 - v_j) < \\varepsilon$:\n",
" - сэмплируется $v_i \\sim \\mathrm{Beta}_{1, \\alpha}(v_i)$,\n",
" - сэпмлируется $\\theta_i \\sim G$,\n",
" - вычисляется соответствующая $\\theta_i$ вероятность $\\pi_i = \\prod_{j=1}^{i-1}(1 - v_j) \\, v_i$,\n",
" - счётчик $i$ увеличивается на 1."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы",
"MCMC"
]
},
"source": [
"## Разделение смеси неопределённого количества распределений\n",
"\n",
"#### Формализация задачи\n",
"\n",
"Пусть имеется смесь неопределённого количества распределений, такая что $i$-е распределение параметризовано вектором $\\theta_i \\in \\mathbb{R}^d$ и имеет вид $p(x \\, \\vert \\, \\theta_i)$, где зависимость плотности от $x$ и $\\theta_i$ дана в аналитическом виде. Обозначим за $\\theta$ бесконечную последовательность неизвестных параметров $\\theta_i$, $\\theta = \\{\\theta_i\\}_{i=1}^{+\\infty}$, а за $\\pi$ обозначим бесконечную последовательность вероятностей соответствующих компонент смеси, $\\pi = \\{\\pi_i\\}_{i=1}^{+\\infty}$. Если на самом деле количество компонент смеси конечное, можно считать, что с какого-то момента все $\\pi_i$ становятся равны 0.\n",
"\n",
"Чисто теоретически, задачу можно бы было поставить как максимизацию по $\\theta$ и $\\pi$ следующего правдоподобия:\n",
"$$p(X, Z \\, \\vert \\, \\theta, \\pi) = \\prod_{i=1}^{l}p(x_i \\, \\vert \\, \\theta_{z_i})p(z_i \\, \\vert \\pi) = \\prod_{i=1}^{l}p(x_i \\, \\vert \\, \\theta_{z_i})\\pi_{z_i},$$\n",
"где $X$ — матрица наблюдений размером $l \\times n$, $l$ — количество объектов, $n$ — количество признаков, $Z$ — вектор размера $l \\times 1$ номеров компонент смеси, породивших соответствующие наблюдения, $x_i$ — признаковое описание $i$-го объекта, $z_i$ — номер компоненты, породившей $i$-й объект. На самом деле, поскольку вектор $Z$ не дан, речь идёт о максимизации не полного правдоподобия, а частичного правдоподобия (правдоподобия только $X$). Однако проблема не в этом. С такой постановкой проблема в том, что она допускает вырожденное решение: можно каждый объект поместить в свою собственную компоненту. Для борьбы с этим часто считают, что количество компонент смеси зафиксировано и что оно меньше количества объектов. Но что делать, если хочется иметь возможность варьировать число компонент смеси по мере поступления новых данных?\n",
"\n",
"Чтобы избежать вырожденных решений, вводится своего рода регуляризция, а именно априорное распределение на пару $(\\theta, \\pi)$. Пара $(\\theta, \\pi)$ может быть проинтерпретирована как вероятностная мера на $\\mathbb{R}^d$. Следовательно, на ум тут же приходит взять в качестве априорного распределения [процесс Дирихле](__home_url__/notes/Распределение Дирихле и случайные процессы Дирихле). Итак, пусть исходя из априорных знаний были подобраны мера $G$ на $\\mathbb{R}^d$ и число $\\alpha > 0$, которыми параметризуется процесс Дирихле $\\xi_{G, \\alpha}$. Будем считать, что именно из такого $\\xi_{G, \\alpha}$ сэмплируется пара $(\\theta, \\pi)$. Тогда может быть выписана уже не условная, а совместная плотность:\n",
"$$p(X, Z, \\theta, \\pi) = \\prod_{i=1}^{l}\\left(p(x_i \\, \\vert \\, \\theta_{z_i})\\pi_{z_i}\\right)p(\\theta, \\pi),$$\n",
"где $p(\\theta, \\pi)$ — вероятность реализации $(\\theta, \\pi)$ случайного процесса Дирихле $\\xi_{G, \\alpha}$. Как это видно из схемы сэмплирования из процесса Дирихле, называемой stick breaking, $p(\\theta, \\pi)$ факторизуется до $p(\\theta)p(\\pi)$, где $p(\\theta) = \\prod_{i=1}^{+\\infty}G(\\theta_i)$.\n",
"\n",
"Задача заключается в том, что по выборке $X$ и выписанному аналитическому виду $p(X, Z, \\theta, \\pi)$ требуется найти апостериорное распределение $p(Z, \\theta, \\pi \\, \\vert \\, X)$.\n",
"\n",
"#### Сколлапсированное сэмплирование Гиббса\n",
"\n",
"Как это часто и бывает в байесовских методах, поставленная задача не решается аналитически. Один из численных способов поиска приближённого решения опирается на марковские цепи Монте-Карло. А поскольку MCMC сходятся долго, хочется сократить количество переменных. Отсюда и возникает коллапсирование. По переменной $\\pi$ проводится схлопывание, и выбрана $\\pi$ просто по той причине, что интеграл возьмётся аналитически:\n",
"$$p(X, Z, \\theta) = \\int p(X, Z, \\theta, \\pi) d\\pi = \\int \\prod_{i=1}^{l}\\left(p(x_i \\, \\vert \\, \\theta_{z_i})p(z_i \\, \\vert \\pi)\\right) p(\\theta)p(\\pi) d\\pi = \\prod_{i=1}^{l}p(x_i \\, \\vert \\theta_{z_i}) \\; p(\\theta) p(Z),$$\n",
"где плотности образующих смесь распределений $p(x_i \\, \\vert \\theta_{z_i})$ были даны исходно, $p(\\theta)$ определена выбором параметра процесса Дирихле $G$, а $p(Z) = \\int \\prod_{i=1}^{l}p(z_i \\, \\vert \\pi)p(\\pi)d\\pi = \\int \\prod_{i=1}^{l}\\pi_{z_i}p(\\pi)d\\pi$.\n",
"\n",
"После коллапсирования сэмплирование Гиббса позволит генерировать пары $(Z, \\theta)$ из распределения $p(Z, \\theta \\, \\vert \\, X)$. При этом сама конструкция сэмплирования Гиббса в этой задаче окажется похожей на схему китайских ресторанов.\n",
"\n",
"Величина $(Z, \\theta)$, которую можно рассмотреть как конкатенацию вектора $Z$ длины $l$ и бесконечной последовательности $\\theta$, формально является бесконечномерной, так что даже один шаг сэмплирования Гиббса должен бы длиться бесконечно долго. Чтобы этого избежать, $\\theta$ понимается как вектор переменной длины, такой что в него входят лишь параметры тех компонент, к которым приписан хотя бы один объект. Если на каком-либо шаге из компоненты перенезначаются в другие компоненты все её объекты, параметры этой компоненты удаляются из вектора $\\theta$.\n",
"\n",
"Итак, пусть схема Гиббса была инициализирована каким-то начальным значением $(Z, \\theta)$. Тогда каждое следующее сэмплирование проводится в три этапа, описываемых ниже.\n",
"\n",
"На первом этапе для каждого из $i$ от 1 до $l$ фиксируются все $z_j$, $j \\ne i$, и весь конечномерный вектор $\\theta$. Обозначим за $Z_{-i}$ вектор длины $l-1$, составленный всеми компонентами $Z$ кроме $i$-й. Нужно просэмплировать новое значение $z_i$ из одномерного распределения $p(z_i \\, \\vert \\, Z_{-i}, \\theta, X)$, известного с точностью до константы. Делается это так:\n",
"* $z_i := k$, где $k$ от 1 до текущей длины вектора $\\theta$, с вероятностью, пропорциональной $n_k^{-i}p(x_i \\, \\vert \\theta_k)$, где $n_k^{-i}$ — количество на текущий момент приписанных к $k$-й компоненте объектов без учёта $i$-го объекта (если он к ней приписан). На то, чему вероятность пропорциональна, можно посмотреть как на популярность компоненты, скорректированную на правдоподобие объекта относительно неё.\n",
"* $z_i := k_*$, где $k_*$ — увеличенная на 1 длина вектора $\\theta$, с вероятностью, пропорциональной $\\alpha\\int p(x_i \\, \\vert \\, \\theta) G(\\theta) d\\theta$. Если выпадает этот исход, длина вектора $\\theta$ увеличивается на 1, и в его конец дописывается компонента $\\theta_{k_*}$, просэмплированная из распределения с плотностью $G(\\theta_{k_*})p(x_i \\, \\vert \\, \\theta_{k_*})$, понимаемого как априорное распределение $G$, перевзвешенное на эмпирические данные.\n",
"\n",
"На втором этапе для каждого из $i$ от 1 до текущей длины $\\theta$ фиксируются все $\\theta_j$, $j \\ne i$, и весь вектор $Z$. Процедура же обновления каждого $\\theta_i$ выглядит так: новое значение $\\theta_i$ сэмплируется из распределения с плотностью $G(\\theta_i) \\prod_{j: z_j = i}p(x_j \\, \\vert \\, \\theta_i)$. \n",
"\n",
"На третьем этапе вектор $\\theta$ модифицируется, а именно из него удаляются все его компоненты, такие что у соответствующих им распределений после первого этапа не осталось приписанных объектов, и оставшиеся компоненты перенумеровываются."
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"байесовские_методы",
"естественные_языки",
"графические_модели",
"тематическое_моделирование"
]
},
"source": [
"## Латентное размещение Дирихле\n",
"\n",
"#### Введение\n",
"\n",
"Латентное размещение Дирихле (Latent Dirichlet Allocation, LDA) является одним из методов тематического моделирования, опирающегося на подход «мешка слов». В рамках такого подхода порядок слов внутри документов не учитывается, а учитывается только то, сколько раз то или иное слово вошло в тот или иной документ.\n",
"\n",
"LDA является довольно простым методом тематического моделирования, и это даёт ему высокую интерпретируемость, так что можно провести такую аналогию, что LDA для задачи тематического моделирования является примерно тем же, чем линейные модели для задачи регрессии.\n",
"\n",
"#### Вероятностная модель\n",
"\n",
"Здесь и далее рассматривается модель без априорного распределения на матрицу тем, а про вариант с априорным распределением будет написано ближе к концу. В рассматриваемом варианте есть следующие сущности:\n",
"* $W$ — наблюдаемый набор документов; у $W$ есть два индекса, но это не матрица, потому что её строки могут иметь разную длину; запись $W_{di}$ пусть обозначает идентификатор слова (далее для краткости будет писаться, что просто слово), стоящего в $d$-м документе на $i$-й позиции; $W_{di} \\in \\{1, \\dots, \\vert V \\vert\\}$, где $V$ — множество всех слов, а запись $\\vert \\cdot \\vert$ применительно ко множеству обозначает количество элементов в нём. Таким образом, чтобы элемент $W_{di}$ существовал, необходимо и достаточно, чтобы индекс $d$ имел значение от 1 до $\\vert D \\vert$, где $D$ — множество всех документов, а индекс $i$ в зависимости от значения индекса $d$ имел значение от 1 до количества слов в $d$-м документе.\n",
"* $Z$ — ненаблюдаемая информация о темах, из которых порождались слова; у $Z$ есть два индекса, и множество допустимых значений этой пары индексов такое же, как у $W$; запись $Z_{di}$ пусть обозначает идентификатор темы (далее для краткости будет писаться, что просто тему), из которой в $d$-м документе было порождено слово, стоящее на $i$-й позиции; $Z_{di} \\in \\{1, \\dots, \\vert T \\vert\\}$, где $T$ — множество всех тем, размер которого $\\vert T\\vert$ является гиперпараметром и должен быть подобран.\n",
"* $\\Theta$ — неизвестная матрица тематических профилей документов, её размеры $\\vert D\\vert \\times \\vert T \\vert$; элемент этой матрицы $\\Theta_{dt}$ показывает вероятность встретить тему $t$ в документе $d$; как следствие, эта матрица не может быть произвольной: на неё накладываются ограничения, что $\\Theta_{dt} \\ge 0$ и для любой $d$-й строки $\\sum_{t=1}^{\\vert T \\vert}\\Theta_{dt} = 1$.\n",
"* $\\Phi$ — неизвестная матрица тем, её размеры $\\vert T \\vert \\times \\vert V\\vert$; элемент этой матрицы $\\Phi_{tw}$ показывает вероятность встретить слово $w$ в теме $t$; как следствие, эта матрица не может быть произвольной: на неё накладываются ограничения, что $\\Phi_{tw} \\ge 0$ и для любой $t$-й строки $\\sum_{w=1}^{\\vert V \\vert}\\Phi_{tw} = 1$.\n",
"* $\\alpha$ — вектор длины $\\vert T \\vert$, такой что все его компоненты строго больше нуля; этот вектор является гиперпараметром и должен быть подобран.\n",
"\n",
"Говоря верхнеуровнево, вероятностная модель для LDA устроена так:\n",
"$$p(W, Z, \\Theta \\, \\vert \\, \\Phi, \\alpha) = \\prod_{d=1}^{\\vert D \\vert}\\left(p(\\Theta_{d\\cdot} \\, \\vert \\, \\alpha) \\prod_{i=1}^{n_d} \\left(p(Z_{di} \\, \\vert \\, \\Theta_{d\\cdot}) p(W_{di}\\, \\vert \\, \\Phi, Z_{di})\\right) \\right),$$\n",
"где $\\Theta_{d\\cdot}$ обозначает $d$-ю строку матрицы $\\Theta$, а $n_{d}$ — количество слов в $d$-м документе. Выписанная формула концептуально интерпретируется так: документы порождаются независимо друг от друга, для каждого документа сначала порождается его тематический профиль, затем в документе поочерёдно порождаются слова, причём для каждого слова сначала порождается его тема, а уже потом из этой темы порождается само слово. Кстати, такую логику можно перерисовать в виде графической модели.\n",
"\n",
"Чтобы специфицировать вероятностную модель для LDA, определим $p(\\Theta_{d\\cdot} \\, \\vert \\, \\alpha)$ как [распределение Дирихле](__home_url__/notes/Распределение Дирихле и случайные процессы Дирихле) с параметром $\\alpha$. Также заметим, что остальные вероятности можно записать в аналитическом виде исходя из определений введённых сущностей:\n",
"$$p(Z_{di} \\, \\vert \\, \\Theta_{d\\cdot}) = \\Theta_{dZ_{di}},$$\n",
"$$p(W_{di}\\, \\vert \\, \\Phi, Z_{di}) = \\Phi_{Z_{di}W_{di}}.$$\n",
"Теперь вероятностную модель для LDA можно записать в явном аналитическом виде:\n",
"$$p(W, Z, \\Theta \\, \\vert \\, \\Phi, \\alpha) = \\prod_{d=1}^{\\vert D \\vert}\\left(\\mathrm{Dir}_{\\alpha}(\\Theta_{d\\cdot}) \\prod_{i=1}^{n_d}\\left(\\Theta_{dZ_{di}} \\Phi_{Z_{di}W_{di}}\\right)\\right) = \\prod_{d=1}^{\\vert D \\vert}\\left(\\frac{\\Gamma(\\sum_{t=1}^{\\vert T \\vert}\\alpha_t)}{\\prod_{t=1}^{\\vert T \\vert}\\Gamma(\\alpha_t)}\\prod_{t=1}^{\\vert T \\vert}\\Theta_{dt}^{\\alpha_t - 1} \\cdot \\prod_{i=1}^{n_d} \\prod_{t=1}^{\\vert T \\vert} \\left(\\Theta_{dt}^{[Z_{di} = t]} \\Phi_{tW_{di}}^{[Z_{di} = t]}\\right)\\right),$$\n",
"где квадратные скобки соответствуют нотации Айверсона.\n",
"\n",
"#### Обучение\n",
"\n",
"При обучении LDA важен не столько сам гиперпараметр $\\alpha$, сколько его длина, то есть количество тем. Если все компоненты $\\alpha$ достаточно малы, то обучение пройдёт успешно. Поэтому будет считать, что гиперпараметр $\\alpha$ сразу подобран удачно и оптимизировать его не нужно.\n",
"\n",
"Обучение LDA сводится к тому, что, имея корпус текстов $W$ и какой-то вектор $\\alpha$, необходимо найти оптимальную структуру тем, то есть матрицу $\\Phi$, количество строк которой задано благодаря $\\alpha$ и количество столбцов которой задано благодаря $W$:\n",
"$$p(W \\, \\vert \\, \\Phi, \\alpha) \\to \\max_{\\Phi}.$$\n",
"Это задача максимизации неполного правдоподобия, ведь правдоподобие ненаблюдаемых переменных $Z$ и $\\Theta$ не учитывается. На самом деле, получится так, что оценки для этих переменных в качестве полезного побочного эффекта будут даны в ходе поиска решения.\n",
"\n",
"Как решать поставленную задачу? Задачи максимизации неполного правдоподобия часто решают при помощи [EM-алгоритма](__home_url__/notes/Общий взгляд на EM-алгоритм). Здесь он тоже позволит получить решение, но с одним нюансом. Точный EM-алгоритм оказывается неприменим, потому что отсутствует сопряжённость по кортежу скрытых переменных, не относящися к параметрам (то есть по объединению $Z$ и $\\Theta$). Из-за этого на Е-шаге задача статистического вывода не решается аналитически. Однако вариационный EM-алгоритм будет работать, потому что условная сопряжённость по кортежу скрытых переменных, не относящихся к параметрам, есть, и, значит, приближённые E-шаги можно делать при помощи mean-field approximation.\n",
"\n",
"Убедимся в этом.\n",
"\n",
"Кортеж наблюдаемых переменных состоит только из $W$, больше наблюдаемых переменных нет. Кортеж скрытых переменных состоит из $\\Phi$, $Z$ и $\\Theta$. Из них $\\Phi$ относится к параметрам, а $Z$ и $\\Theta$ нет. Значит, на E-шаге при текущем зафиксированном значении $\\overline{\\Phi}$ будет оцениваться апостериорное распределение $p(Z, \\Theta \\, \\vert \\, W, \\overline{\\Phi}, \\alpha)$. В вариационном EM-алгоритме оценка ищется в виде $q(Z)q(\\Theta)$, и искусственно вводимая таким выбором факторизация приводит к приближённости.\n",
"\n",
"Теперь перейдём непосредственно к проверке условной сопряжённости по кортежу скрытых переменных, не относящихся к параметрам (в [заметке про mean-field approximation](__home_url__/notes/Вариационный байесовский вывод и mean-field approximation) она называлась условная $T$-сопряжённость, потому что в той заметке символ $T$ обозначал не множество тем, а кортеж скрытых переменных, не относящихся к параметрам). Применительно к LDA определение сводится к наличию следующих свойств:\n",
"* При зафиксированных значениях $\\overline{\\Phi}$ и $\\overline{\\Theta}$ условное априорное распределение $p(Z \\, \\vert \\, \\overline{\\Theta}, \\overline{\\Phi}, \\alpha)$ является сопряжённым к условному правдоподобию $p(W \\, \\vert \\, Z, \\overline{\\Theta}, \\overline{\\Phi}, \\alpha)$.\n",
"* При зафиксированных значениях $\\overline{\\Phi}$ и $\\overline{Z}$ условное априорное распределение $p(\\Theta \\, \\vert \\, \\overline{Z}, \\overline{\\Phi}, \\alpha)$ является сопряжённым к условному правдоподобию $p(W \\, \\vert \\, \\Theta, \\overline{Z}, \\overline{\\Phi}, \\alpha)$.\n",
"\n",
"Для примера проверим первое свойство, второе проверяется аналогично. Условное априорное распределение $p(Z \\, \\vert \\, \\overline{\\Theta}, \\overline{\\Phi}, \\alpha)$ имеет плотность, принадлежащую следующему классу функций от $Z$: $\\prod_{d=1}^{\\vert D \\vert} \\prod_{i=1}^{n_d} \\prod_{t=1}^{\\vert T \\vert} \\mathrm{const}^{[Z_{di} = \\mathrm{const}]}$. После умножения этой плотности на условное правдоподобие $p(W \\, \\vert \\, Z, \\overline{\\Theta}, \\overline{\\Phi}, \\alpha) = \\prod_{d=1}^{\\vert D \\vert} \\prod_{i=1}^{n_d} \\prod_{t=1}^{\\vert T \\vert} \\left(\\overline{\\Phi}_{tW_{di}}^{[Z_{di} = t]}\\right)$ произведение останется в том же самом классе функций от $Z$. Значит, сопряжённость имеется.\n",
"\n",
"Тем самым было показано, что при применении вариационного EM-алгоритма к задаче обучения LDA никаких непреодолимых вычислительных препятствий не возникнет, и решение можно будет найти.\n",
"\n",
"#### Модификации\n",
"\n",
"В вероятностную модель для LDA можно включить априорное распределение $\\Phi$. Для этого добавляют ещё один гиперпараметр: $\\beta$, вектор длины $\\vert V \\vert$, которым параметризуется распределение Дирихле, из которого порождаются строки матрицы $\\Phi$.\n",
"\n",
"Вероятностная модель в таком варианте принимает вид:\n",
"$$p(W, \\Phi, Z, \\Theta \\, \\vert \\, \\alpha, \\beta) = \\left(\\prod_{t=1}^{\\vert T \\vert} p(\\Phi_{t\\cdot} \\, \\vert \\, \\beta)\\right)\\prod_{d=1}^{\\vert D \\vert}\\left(p(\\Theta_{d\\cdot} \\, \\vert \\, \\alpha) \\prod_{i=1}^{n_d} \\left(p(Z_{di} \\, \\vert \\, \\Theta_{d\\cdot}) p(W_{di}\\, \\vert \\, \\Phi, Z_{di})\\right) \\right),$$\n",
"где $p(\\Phi_{t\\cdot} \\, \\vert \\, \\beta) = \\mathrm{Dir}_{\\beta}(\\Phi_{t\\cdot})$.\n",
"\n",
"Отдельно стоит заметить, что обучать LDA можно не только вариационным EM-алгоритмом, но и при помощи [сколлапсированного сэмплирования Гиббса](__home_url__/notes/Разделение смеси неопределённого количества распределений)."
]
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 2
}