Федеральное государственное автономное образовательное учреждение
высшего профессионального образования
«СИБИРСКИЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ»

arXiv:1601.07224v1 [] 26 Jan 2016

Институт математики и фундаментальной информатики
Базовая кафедра математического моделирования и процессов управления

БАКАЛАВРСКАЯ РАБОТА
Направление 010101.62 Математика

ПОРОЖДАЮЩЕЕ ВЕРОЯТНОСТНОЕ ПРОГРАММИРОВАНИЕ

Выпускник – Юрий Перов1.
Научный руководитель в СФУ – Татьяна Валерьевна Крупкина.

Кембридж-Оксфорд-Красноярск 2012–2014

1

Адрес для корреспонденции: yuraperov@gmail.com

РЕФЕРАТ
Выпускная квалификационная бакалаврская работа по теме «Порождающее вероятностное программирование» содержит 48 страниц текста, 50 использованных источников, 12 рисунков.
Ключевые слова: ВЕРОЯТНОСТНОЕ ПРОГРАММИРОВАНИЕ, МАШИННОЕ
ОБУЧЕНИЕ, ИСКУССТВЕННЫЙ ИНТЕЛЛЕКТ, ВЕРОЯТНОСТНЫЕ МОДЕЛИ, АВТОМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ.
Работа посвящена новому направлению в области машинного обучения и компьютерных наук — вероятностному программированию. В работе дается краткое реферативное введение в языки вероятностного программирования Church/Venture/Anglican, а также описываются результаты первых экспериментов по автоматической генерации вероятностных программ.

СОДЕРЖАНИЕ
Введение

4

1 Краткое введение в языки вероятностного программирования Church,
Venture и Anglican

8

1.1

Первое знакомство с Church, Venture, Anglican . . . . . . . . . . . . . . . . .

8

1.2

Статистический вывод в вероятностных языках программирования с помощью алгоритма Метрополиса-Гастингса . . . . . . . . . . . . . . . . . . . . . 11
1.2.1

Метод «выборки с отклонением» . . . . . . . . . . . . . . . . . . . . . 11

1.2.2

Пространство историй выполнений вероятностных программ . . . . . 12

1.2.3

Апостериорное распределение историй выполнений программ . . . . . 14

1.2.4

Использование методов Монте-Карло по схеме Марковских цепей . . 15

1.2.5

Программная реализация статистического вывода . . . . . . . . . . . 19

1.3

Эффективность вывода . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

1.4

Порождающее вероятностное программирование в распознавании образов . 24

1.5

О различиях между Church, Venture, Anglican . . . . . . . . . . . . . . . . . . 26

2 Автоматическая генерация вероятностных программ

27

2.1

Обзор литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.2

Описание подхода . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.3

Грамматика и порождающие правила . . . . . . . . . . . . . . . . . . . . . . . 33

2.4

Вероятности использования порождающих правил . . . . . . . . . . . . . . . 35

2.5

Эксперименты . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.1

Выборки из сгенерированных вероятностных программ . . . . . . . . 36

2.5.2

Вывод вероятностных программ, определяющих распределения, совпадающие или приближающие классические одномерные распределения . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2

2.5.3

Обобщение произвольных эмпирических данных с помощью порождающего вероятностного программирования . . . . . . . . . . . . . . . 37

2.5.4
2.6

«Компиляция» вероятностных программ . . . . . . . . . . . . . . . . . 39

Обобщение порождающего вероятностного программирования . . . . . . . . 41

Заключение

44

Список использованных источников

45

3

Моей любимой Маше и всем
другим моим любимым родным

ВВЕДЕНИЕ
Данная бакалаврская работа посвящена вероятностному программированию [1], новому направлению в областях машинного обучения, искусственного интеллекта и компьютерных наук, а именно реферативному краткому введению в вероятностное программирование, описанию языков вероятностного программирования Church [2], Venture [3]
и Anglican [4], описанию подхода порождающего вероятностного программирования для
решения задач распознавания образов [5], а также представлению полученных предварительных результатов по автоматизации вывода вероятностных моделей для вероятностного программирования [6].
Описание структуры работы
Первая часть данной работы начинается с общих сведений о вероятностном программировании, а затем в ней кратко описываются результаты, полученные коллегами автора и им
самим в работе над научными проектами на протяжении двух лет в Массачусетском технологическом инситуте под руководством профессора Джошуа Тененбаума и доктора Викаша Мансингхи (Кембридж, штат Массачусетс, США) и Оксфордском университете под
руководством профессора Френка Вуда (Великобритания, Оксфорд). Эта часть работы,
являясь по сути переводом, представляет собой реферативную выдержку о вероятностном программировании и его приложениях на русском языке. Насколько автору известно,
литературы о вероятностном программировании на русском языке практически нет, поэтому он надеется, что данная чисто реферативная часть работы принесет существенную
пользу русскоговорящему научному сообществу, а особенно заинтересованным студентам
и школьникам, которые впервые захотят познакомиться с развивающимся направлением
вероятностного программирования.
Во второй части данной работы автор описывает новые результаты в области вероятностного программирования, связанные с автоматизированным или полуавтоматизированным выводом вероятностных моделей для вероятностного программирования, полученные во время стажировки автора в Департаменте технических наук Оксфордского

4

университета в научной лаборатории профессора Френка Вуда и под его руководством.
Вторая часть завершается рассуждениями автора об обобщении автоматизированного изучения и вывода компьютером вероятностных моделей в виде вероятностных программ, то
есть о возможностях порождающего вероятностного программирования.
Вероятностное программирование
Вероятностное программирование можно определить как компактный, композиционный
способ представления порождающих вероятностных моделей и проведения статистического вывода в них с учетом данных с помощью обобщенных алгоритмов.
Вероятностная модель является важным понятием машинного обучения, одного из
основных направлений искусственного интеллекта на сегодняшний день. В общем случае в рамках теории машинного обучения перед компьютером ставится задача произвести
какое-то действие yi на основе входных данных xi , априорных знаний и возможности взаимодействовать со средой. Без ограничений общности целевым действием для компьютера
можно считать производство ответа в виде выходных данных, представленных в виде информации. Например, в робототехнике эта информация может являться инструкциями
для моторов и механических устройств для выполнения тех или иных физических действий роботом.
При использовании различных подходов в рамках машинного обучения используются модели, которые являются формальным «описанием компонентов и функций, отображающих существенные свойства моделируемого объема или процесса» [7]. В рамках
машинного обучения используются вероятностные модели, так как свойства, элементы и
связи между ними являются не фиксированными, а стохастическими.
Классическим примером одного из подходов в машинном обучении является способ
«обучение с учителем», когда существует N прецедентов, то есть пар входных и выходных
данных обучающей выборки X = {xi }, Y = {yi }, i = 1, N , и необходимо найти алгоритм
описывающий зависимость между xi и yi , то есть алгоритм F , который позволяет на основе каждого элемента входных данных xi получать абсолютно или достаточно точный
элемент выходных данных yi , то есть F (xi ) → yi . С помощью данного алгоритма затем
для M известных элементов входных данных {xi }, j = N + 1, N + M находятся значения
M неизвестных элементов выходных данных {yi }, j = N + 1, N + M . В рамках машинного обучения эта проблема, в том числе, решается с помощью методов регрессионного
5

анализа.
В только что приведенной как пример модели неизвестными (скрытыми) параметрами T = {tj }, j = 1, K будет информация о характеристиках модели F , которые подлежат
выводу при данной обучающей выборке в виде пар (xi , yi ).
Приведем простой пример: в линейной регрессии (при L независимых переменных)
значениями независимых переменных будут xi,j , значениями зависимой переменной будут
yi , параметрами модели будут t1 , t2 , . . . , tL+1 , а алгоритмом F будет
yi = Ft1 ,...,tL+1 (xi ) = t1 + t2 · xi,1 + . . . + tL+1 · xi,L .
Для простоты без ограничений общности будем считать, что мы можем всегда получить yi детерминированно, зная T и xi .
В порождающих вероятностных моделях задается совместное распределение вероятностей P (T, X), обычно сначала путем задания априорного распределения P (T ), а затем
задания условного распределения P (X | T ). Это и называется моделью.
При заданной модели и известных X задачей будет являться поиск апостериорного распределения на T , таким образом P (T | X). Одним из способов поиска данного
апостериорного распределения является применение теоремы Байеса:

P (T | X) =
где P (X) теоретически можно найти как

R

P (T )P (X | T )
,
P (X)
P (T )P (X | T ) dT , но при решении практиче-

ских задач это часто невозможно, так как перебор всего пространства T не поддается
аналитическому решению или решению с помощью численных методов за разумное время. Поэтому чаще всего при решении задач в рамках Байесовского подхода стоит задача
поиска ненормированного значения P (T | X):
P (T | X) ∝ P (T )P (X | T ),
где символ ∝, часто встречающийся в зарубежной литературе, но очень редко встречающийся у нас, означает «пропорционально». Нормировочную константу затем можно найти
приближенно, но иногда ее значение даже не вычисляют, так как бывает достаточно найти
и работать дальше с наиболее вероятными элементами T̂ из апостериорного распределения
6

P (T | X).
Обычно, особенно при решении практических задач с большим объемом данных
и в рамках сложных моделей, апостериорное распределение P (T | X) находят не точно,
а с помощью приближенных методов, в том числе с помощью методов Монте-Карло [8],
которые позволяют сгенерировать выборку из интересующего нас распределения.
Как отмечалось в самом начале данного подраздела, вероятностное программирование позволяет:
1. Композиционно и компактно записывать порождающую вероятностную модель
с помощью задания априорных вероятностей P (T ) и P (X | T ) в виде алгоритма
(вероятностной программы) с использованием стохастических функций (например, стохастическая функция Normal(µ, σ)) и вспомогательных детерминированных функций (например, + или ∗).
2. Снабжать модель данными X̂, таким образом теоретически определяя условное
распределение P (T | X̂).
3. Производить статистический вывод для генерации выборки из условного распределения с помощью обобщенных алгоритмов статистического вывода, подходящих для всех или большого множества моделей.
Для ознакомления с машинным обучением и искусственным интеллектом автор
рекомендует следующие источники: [9, 10, 11]. Информацию о вероятностных моделях и
Байесовских методах на русском языке можно найти в [12].

7

1

Краткое введение в языки вероятностного программирования
Church, Venture и Anglican

Существует более 15 языков вероятностного программирования, перечень с кратким описанием каждого из них можно найти на [13]. В данной работе реферативно будут рассмотрены три языка вероятностного программирования: Church [2], Venture [3] и Anglican [4].
Языки Venture и Anglican являются продолжениями языка Church. Church в свою очередь
основан на языке «обычного» программирования Lisp и Scheme. Заинтересованному читателю крайне рекомендуется ознакомиться с книгой [14], являющейся одним из лучших
способов начать знакомство с языком «обычного» программирования Scheme.
1.1

Первое знакомство с Church, Venture, Anglican

На текущий момент любой язык вероятностного программирования напрямую связан с
методами статистического вывода, которые используются в нем, поэтому часто язык ассоциируется с платформой вероятностного программирования, то есть с его технической
реализацией в виде компьютерной программы.
Рассмотрим задание простой вероятностной модели Байесовской линейной регрессии [10] на языке вероятностного программирования Venture/Anglican [15] в виде вероятностной программы:
1 [ASSUME t 1 ( normal 0 1 ) ]
2 [ASSUME t 2 ( normal 0 1 ) ]
3 [ASSUME n o i s e 0 . 0 1 ]
4 [ASSUME noisy_x ( lambda ( time ) ( normal (+ t 1 (∗ t 2 time ) ) n o i s e ) ) ]
5 [OBSERVE ( noisy_x 1 . 0 ) 1 0 . 3 ]
6 [OBSERVE ( noisy_x 2 . 0 ) 1 1 . 1 ]
7 [OBSERVE ( noisy_x 3 . 0 ) 1 1 . 9 ]
8 [ PREDICT t 1 ]
9 [ PREDICT t 2 ]
10 [ PREDICT ( noisy_x 4 . 0 ) ]
Скрытые искомые параметры — значения коэффициентов t1 и t2 линейной функции
x (time) = t1 + t2 · time. У нас есть априорные предположения о данных коэффициентах,
8

а именно мы предполагаем, что они распределены по закону нормального распределения
Normal(0, 1) со средним 0 и стандартным отклонением 1. Таким образом, мы определили в
первых двух строках вероятностной программы вероятность P (T ), описанную в предыдущем раздел. Инструкцию [ASSUME name expression] можно рассматривать как определение случайной величины с именем name, принимающей значение вычисляемого выражение
(программного кода) expression, которое содержит в себе неопределенность.
Вероятностные языки программирования (здесь и далее будут иметься в виду конкретно Church, Venture, Anglican, если не указано иное), как и Lisp/Scheme, являются
функциональными языками программирования, и используют польскую нотацию1 при
записи выражений для вычисления. Это означает, что в выражении вызова функции сначала располагается оператор, а уже только потом аргументы: (+ 1 2), и вызов функции
обрамляется круглыми скобками. На других языках программирования, таких как C++
или Python, это будет эквивалентно коду 1 + 2.
В вероятностных языках программирования выражение вызова функции принято
разделять на три разных вида:
1. Вызов детерминированных процедур (primitive-procedure arg1 ...argN), которые при одних и тех же аргументах всегда возвращают одно и то же значение.
К таким процедурам, например, относятся арифметические операции.
2. Вызов вероятностных (стохастических) процедур (stochastic-procedure arg1
...argN), которые при каждом вызове генерируют случайным образом элемент
из соответствующего распределения. Такой вызов определяет новую случайную
величину. Например, вызов вероятностной процедуры (normal 1 10) определяет случайную величину, распределенную по закону нормального распределения
√
Normal(1, 10), и результатом выполнения каждый раз будет какое-то вещественное число.
3. Вызов

составных

процедур

(compound-procedure arg1 ...argN),

где

compound-procedure — введенная пользователем процедура с помощью специального выражения lambda: (lambda (arg1 ...argN) body), где body — тело
процедуры, состоящее из выражений. В общем случае составная процедура
является стохастической (недетерминированной) составной процедурой, так как
1 Venture имеет отдельный дополнительный вид синтаксиса VentureScript, использующий инфиксную нотацию и не требующий обрамления вызова функций скобками, то есть схожий по своей сути с привычными большинству людей языками
программирования C, C++, Python и т.д.

9

ее тело может содержать вызовы вероятностных процедур.
После этого мы хотим задать условную вероятность P (X | T ) наблюдаемых переменных x1 , x2 , x3 при заданных значениях скрытых переменных t1 , t2 и параметра
time.
Перед вводом непосредственно самих наблюдений с помощью выражения [OBSERVE
...] мы определяем общий закон для наблюдаемых переменных {xi } в рамках нашей
модели, а именно мы предполагаем, что данные наблюдаемые случайные величины при
заданных t1 , t2 и заданном уровне шума noise распределены по закону нормального рас√
пределения Normal(t1 +t2 ·time, noise) со средним t1 +t2 ·time и стандартным отклонением
noise. Данная условная вероятность определена на строках 3 и 4 данной вероятностной
программы. noisy_x определена как функция, принимающая параметр time и возвращающая случайное значение, определенное с помощью вычисления выражение (normal (+
t1 (* t2 time)) noise) и обусловленное значениями случайных величин t1 и t2 и переменной noise. Отметим, что выражение (normal (+ t1 (* t2 time)) noise) содержит
в себе неопределенность, поэтому каждый раз при его вычислении мы будем получать в
общем случае разное значение.
На строках 5—7 мы непосредственно вводим известные значения xˆ1 = 10.3, xˆ2 =
11.1, xˆ3 = 11.9. Инструкция вида [OBSERVE expression value] фиксирует наблюдение о
том, что случайная величина, принимающая значение согласно выполнению выражения
expression, приняла значение value.
Повторим на данном этапе всё, что мы сделали. На строках 1—4 с помощью инструкций вида [ASSUME ...] мы задали непосредственно саму вероятностную модель:
P (T ) и P (X | T ). На строках 5—7 мы непосредственно задали известные нам значения
наблюдаемых случайных величин X с помощью инструкций вида [OBSERVE ...].
На строках 8—9 мы запрашиваем у системы вероятностного программирования апостериорное распределение P (T | X) скрытых случайных величин t1 и t2 . Как уже было
сказано, при большом объеме данных и достаточно сложных моделях получить точное
аналитическое представление невозможно, поэтому инструкции вида [PREDICT ...] генерируют выборку значений случайных величин из апостериорного распределения P (T | X)
или его приближения. Инструкция вида [PREDICT expression] в общем случае генерирует один элемент выборки из значений случайной величины, принимающей значение

10

согласно выполнению выражения expression. Если перед инструкциями вида [PREDICT
...] расположены инструкции вида [OBSERVE ...], то выборка будет из апостериорного
распределения2 , обусловленного перечисленными ранее введенными наблюдениями.
Отметим, что в завершении мы можем также предсказать значение функции
x(time) в другой точке, например, при time = 4.0. Под предсказанием в данном случае понимается генерация выборки из апостериорного распределения новой случайной
величины при значениях скрытых случайных величин t1 , t2 и параметре time = 4.0.
Для генерации выборки из апостериорного распределения P (T | X) в языке
программирования Church в качестве основного используется алгоритм МетрополисаГастингса, который относится к методам Монте-Карло по схеме Марковских цепей. В
следующем подразделе будет произведено подробное описание применения данного алгоритма для обобщенного статистического вывода в вероятностных языках. Под «обобщенным» выводом в данном случае понимается то, что алгоритм может быть применен
к любым вероятностным программам, написанным на данном вероятностном языке программирования.
1.2

Статистический вывод в вероятностных языках программирования с помощью алгоритма Метрополиса-Гастингса

Описание алгоритма Метрополиса-Гастингса в применении к «семейству» вероятностных
языков Church впервые опубликовано в [2] и более подробно описано в [16].
Получить выборку из N элементов из априорного распределения скрытых параметров P (T ), наблюдаемых величин P (X | T ) или их совместного распределение P (T, X)
какой-либо порождающей вероятностной модели, записанной в виде вероятностной программы, не составляет труда. Для этого достаточно выполнить вероятностную программу
N раз. Отметим очевидный факт, что в данном случае вероятностная программа будет содержать лишь инструкции вида [ASSUME ...] (задание вероятностной модели) и [PREDICT
...] (перечисление случайных величин, выборку которых мы генерируем).
1.2.1

Метод «выборки с отклонением»

При заданных наблюдениях с помощью инструкций вида [OBSERVE expression value]
наиболее простым способом получения апостериорного распределения является метод
2 Говоря

точнее, конечно, из приближения апостериорного распределения.

11

«выборки с отклонением» [17]. Для понимания рассмотрим следующую вероятностную
программу:
1 [ASSUME a ( uniform−d i s c r e t e 1 6 ) ]
2 [ASSUME b ( uniform−d i s c r e t e 1 6 ) ]
3 [OBSERVE (+ a b ) 5 ]
4 [ PREDICT (∗ a b ) ]
Отметим, что стохастическая процедура (uniform-continuous a b) возвращает
значение случайной величины, распределенной по равномерному дискретному закону распределения с носителем {a, . . . , b}.
Текстом задачу, записанную выше с помощью вероятностной программы, можно
сформулировать следующим образом: подбрасываются два шестигранных игральных кубика, в сумме выпало 5; каково распределение произведения очков на этих двух кубиках?
Метод «выборки с отклонением» заключается в том, чтобы генерировать значения
очков на первом и втором кубиках из их априорного распределения (таким образом, два
независимых дискретных равномерных распределения) и проверять, равна ли их сумма
5. Если нет, данная попытка отвергается и не учитывается. Если да, то произведение
очков на кубиках добавляется во множество элементов выборки. Данная выборка и будет аппроксимацией апостериорного распределения значений произведения очков на двух
кубиках, если известно, что сумма очков равна 5.
Данный метод является неэффективным, и при решении более сложных задач просто вычислительно неосуществимым. Также отметим, что он хорошо подходит только для
дискретных значений случайных и промежуточных переменных.
1.2.2

Пространство историй выполнений вероятностных программ

При выполнении вероятностной программы каждая случайная величина принимает определенное значение. Например, следующая вероятностная программа содержит три случайных величины, каждая из которых распределена по закону нормального распределения:
1 [ASSUME a ( normal 0 1 ) ]
2 [ASSUME b ( normal 0 1 ) ]

12

3 [ASSUME c ( normal (+ a b ) 1 ) ]
4 [ PREDICT c ]
При выполнении данной вероятностной программы каждая из трех случайных величин примет свое случайное значение. Назовем историей выполнения вероятностной
программы отображение, сопоставляющее каждой случайной величине ее значение. Мощность истории выполнения совпадает с количеством случайных величин, которые были
созданы при ее выполнении.
Для каждой программы можно определить соответствующее вероятностное пространство, элементарными событиями которого будут являться все истории ее выполнения. Вероятностью каждого элементарного события будет являться произведение вероятностей принятие того или иного значения каждой случайные величины при данной
истории выполнения.
Очевидно, что история выполнения однозначным образом определяет выполнение
вероятностной программы и принимаемые значения, как случайных, так и зависящих от
них детерминированных выражений.
Для примера рассмотрим еще более простую вероятностную программу.
1 [ASSUME a ( b e r n o u l l i 0 . 7 ) ]
2 [ASSUME b ( b e r n o u l l i 0 . 7 ) ]
3 [ PREDICT a ]
4 [ PREDICT b ]
В ней обе случайных величины имеют свое название, что облегчает запись историй
выполнений (иначе необходимо вводить адресную схему для наименования случайных
величин, чтобы однозначно идентифицировать их).
Данная вероятностная программа имеет четыре различных возможных историй
выполнений, а именно {a : 0, b : 0}, {a : 0, b : 1}, {a : 1, b : 0} и {a : 1, b : 1} с вероятностями
0.09, 0.21, 0.21 и 0.49 соответственно.
Отметим также, что количество «активных» случайных величин при выполнении
одной и той же вероятностной программы может меняться, как, например, в следующей
программе:
[ASSUME g e o m e t r i c

13

( lambda ( p ) ( i f (= ( b e r n o u l l i p ) 1 ) 0 (+ ( g e o m e t r i c p ) 1 ) ) ) ]
[ PREDICT ( g e o m e t r i c 0 . 5 ) ]
Данная программа генерирует элемент (т.е. одноэлементную выборку) из геометрического распределения с параметром 0.5 с помощью составной процедуры geometric,
параметризованной параметром p. При выполнении данной вероятностной программы
в истории ее выполнения количество «активных» случайных величин может быть любым.
1.2.3

Апостериорное распределение историй выполнений программ

При добавлении наблюдений, то есть фиксации значения определенных случайных величин, можно считать, что рассматривается подмножество множества элементарных событий, а именно те элементы, в историях выполнений которых наблюдаемые случайные
величины принимают желаемое значение.
Можно описать другое вероятностное пространство, множеством элементарных событий которого будет являться только что описанное подмножество. Вероятностная мера
в новом вероятностном пространстве может быть «индуцирована» вероятностной мерой из
первоначального вероятностного пространства с учетом нормировочной константы.
Только что описанный переход полностью сочетается с теоремой Байеса, которая в
нашем случае записывается следующим образом:

P (T | X = x̂) =

P (T )P (X = x̂ | T )
,
P (X = x̂)

где X — множество случайных величин, для которых мы знаем фиксированные наблюдаемые значения; T — множество случайных величин, ассоциируемых со скрытыми параметрами, апостериорное распределение которых мы заинтересованы получить (грубо
говоря, это всё остальные случайные величины); x̂ — наблюдаемые значения случайных
величин.
Отметим, что часто мы заинтересованы не в апостериорном распределении
P (T | X = x̂), а в апостериорном распределении лишь части скрытых параметров
P (T 0 ⊂ T | X = x̂) или даже функции f от них P (f (T ) | X = x̂), в общем случае не
являющейся биекцией. С другой стороны, так как мы не ставим задачу получить аналитическое представление данных апостериорных распределений, а лишь выборку из них, то
14

это не играет большой роли в нашем случае при решении задач методами Монте-Карло:
мы можем генерировать элементы выборки из P (T | X = x̂), и затем использовать только
значения нужных нам скрытых случайных величин и/или действовать функцией f на
них.
1.2.4

Использование методов Монте-Карло по схеме Марковских цепей

Математическо-статистический аппарат методов Монте-Карло по схеме Марковских цепей
кратко и «современно» изложен в [8].
Интуитивно опишем, что мы собираемся делать. У нас есть вероятностная программа g, определяющая множество скрытых T и наблюдаемых X случайных величин.
Вероятностная программа своей записью задает априорное распределение P (T ), P (X | T ),
а значит и совместное распределение P (T, X). Каждый раз путем выполнения данной
программы мы получаем одну из возможных реализаций данной программы, которая биективно описывается соответствующей историей выполнения hi ∈ H, где множество H —
множество всех возможных историй выполнений данной вероятностной программы, которое можно рассматривать как множество элементарных событий. Для каждой истории
выполнения определена ее вероятностная мера
P (hi ) = P (T = t̃, X = x̃) = P (T = t̃)P (X = x̃ | T = t̃).
У нас также есть «экспериментальное» значение каждой наблюдаемой случайной
величины x̂i , то есть множество {x̂i } = x̂. Существует подмножество историй выполнений
H 0 ⊂ H, в котором наблюдаемые случайные величины принимают желаемое значение.
Данное подмножество можно рассматривать как множество элементарных событий другого вероятностного пространства, вероятностную меру на котором можно «индуцировать»
из предыдущего путем деления на нормирующую постоянную P (X = x̂):
P (h0 ) =

P (T = t̃)P (X = x̂ | T = t̃)
.
P (X = x̂)

Будем обозначать ненормированную вероятностную меру
˜ 0 ) = P (T = t̃)P (X = x̂ | T = t̃).
P (h

15

Рассмотрим цепь Маркова, исходами которой являются истории выполнений h0 ∈
H 0 . В качестве начального распределения можно выбрать априорное распределение
P (T )P (X | T ), при котором значения наблюдаемых случайных величин не выбираются
случайным образом согласно их распределению, а устанавливаются согласно их значениям. Например, в вероятностной программе
1 [ASSUME a (gamma 1 1 ) ]
2 [ASSUME b ( lambda ( ) ( normal a 1 ) ) ]
3 [OBSERVE ( b ) 5 . 3 ]
первая случайная величина a, распределенная по закону Гамма-распределения
Gamma(1, 1), будет являться скрытой, и будет сгенерирована согласно данному закону
распределения. Наблюдаемая же случайная величина, распределенная по нормальному
закону Normal(a, 1), не будет сгенерирована, а будет установлена в соответствии с ее наблюдаемым значением.
Мы хотим установить такие правила перехода по схеме Метрополиса-Гастингса из
одного состояния (исхода) цепи Маркова в другое, чтобы стационарное распределение
данной цепи Маркова совпадало с распределением P (h0 ). В таком случае для получения
аппроксимации искомого апостериорного распределения в виде выборки нам будет достаточно имитировать данную цепь Маркова [8, 18].
В алгоритме Метрополиса-Гастингса вероятностная мера может быть известна с
точностью до нормировочной константы, что и происходит в нашем случае. На каждом
шаге алгоритма дано текущее состояние h0t и в соответствии с заданным заранее условным
распределение предлагается новое состояние h0∗ ∼ Q( · | h0t ). Таким образом, Q можно назвать распределением предлагаемых переходов. После этого подсчитывается коэффициент
«принятия» нового состояния:
P̃ (h0∗ ) Q(h0t | h0∗ )
α = min 1,
P̃ (h0t ) Q(h0∗ | h0t )

!
.

Состояние h0∗ принимается в качестве следующего состояния h0t+1 с вероятностью
α, в противном случае h0t+1 := h0t .
Новое состояние предлагается следующим образом:
1. Случайным образом (равномерно) выбирается одна «активная» скрытая случайная величина r ∈ dom h0t для «вариации».
16

2. Предлагается новое значение данной случайной величины r∗ ∼ κ( · | r = r̂, h0t ).
Условная вероятность κ в данном случае будет локальным для r распределением
предлагаемых переходов.
3. Если случайная величина r влияет на поток выполнения вероятностной программы, и в результате ее нового значения должны быть исполнены другие ветви
выполнения вероятностной программы, они исполняются и в общем случае происходит генерация новых случайных величин, которые прежде были неактивны.
Для примера рассмотрим следующую вероятностную программу:
1 [ASSUME a ( b e r n o u l l i 0 . 3 ) ]
2 [ASSUME b ( i f (= a 1 ) ( normal 0 1 ) (gamma 1 1 ) ) ]
3 [ PREDICT c ]
В данной вероятностной программе три случайных величины, первая ξ1 распределена по
закону Бернулли, вторая ξ2 по закону нормального распределения, третья ξ3 по закону
Гамма-распределения. ξ1 можно также называть «a» (по имени переменной), хотя переменная «b» не является сама по себе случайной величиной, а зависит от значения ξ1 ,
определяющей поток выполнения вероятностной программы, и от значения либо ξ2 , либо
ξ3 . Отметим также, что при каждом выполнении данной вероятностной программы будет
существовать только две активных случайных величины.
Предположим, что текущим состоянием вероятностной программы в момент времени t была история выполнения h0t = {ξ1 = 1, ξ2 = 3.2}. Это означает, что случайная
величина ξ приняла значение 1, и поэтому для задания переменной b был сгенерирован
элемент из нормального распределения, то есть была «реализована» случайная величина
ξ2 ∼ Normal(0, 1). Очевидно, что
P (h0t ) = P (ξ1 = 1) · P (ξ2 = 3.2 | ξ1 = 1) = 0.3 · fNormal(0,1) (3.2) ≈ 0.0007152264603.
Предложим новое состояние h0∗ :
1. Выберем случайным образом одну из двух «активных» случайных величин,
пусть это будет ξ1 .
2. Предложим случайным образом новое значение данной случайной величины согласно ее априорному распределению (то есть Bernoulli(0.3)), пусть это будет 0.
3. Так как случайная величина ξ1 действительно влияет на поток выполнения ве17

роятностной программы и при изменении ее значения в данном случае должна
быть выполнена другая ветвь, выполним данную ветвь, генерируя значения «активирующихся» случайных величин (в нашем случае только ξ3 ). Предположим,
что ξ3 стало равным 12.3.
Тогда
Q(h0∗ | h0t ) =

1
· 0.7 · fGamma(1,1) (12.3).
2

В общем случае
Q(h0∗ | h0t ) =

1
· κ(r∗ | r = r̂, h0t ) · Rnew ,
0
|ht |

где |h0t | — количество активных случайных величин в текущей истории выполнения,
κ(r∗ | r = r̂, h0t ) — вероятность нового значения варьируемой случайной величины r, Rnew
— совместная вероятность выбора своих значений «активирующихся» случайных величин.
При уже фиксированном h0∗ обратное Q(h0t | h0∗ ) для коэффициента принятия алгоритма Метрополиса-Гастингса в общем случае может быть посчитано аналогичным образом:
Q(h0t | h0∗ ) =

1
· κ(r̂ | r = r∗ , h0∗ ) · Rold ,
|h0∗ |

где Rold — совместная вероятность выбора своих значений случайных величин, активных
в h0t , но не являющихся активными в h0∗ .
В только что рассмотренном примере
Q(h0t | h0∗ ) =

1
· 0.3 · fNormal(0,1) (3.2).
2

И тогда с учетом того, что P (h0∗ ) = 0.7 · fGamma(1,1) (12.3) ≈ 0.000003186221124, мы
получаем, что


0.7 · fGamma(1,1) (12.3) · 12 · 0.3 · fNormal(0,1) (3.2)
α = min 1,
= 1.
0.3 · fNormal(0,1) (3.2) · 21 · 0.7 · fGamma(1,1) (12.3)
Поэтому в данном конкретном примере h0t+1 = h0∗ в любом случае, то есть с вероятностью один.
Выбор локального для r распределения предлагаемых переходов может быть раз-

18

ным. В простейшем случае, если r — независимая случайная величина, данное распределение выбирается идентичным априорному распределению для r. Если же r — перестановочная случайная величина [19], а такие случайные величины поддерживаются рассматриваемыми языками вероятностного программирования, то ее распределение может быть
выбрано с учетом уже накопленных значений.
Таким образом, мы описали алгоритм для предложения нового состояния h0∗ при
текущем состоянии h0t и описали его в рамках метода Метрополиса-Гастингса. Чтобы получить выборку {h0j }N
j=1 из N элементов желаемого апостериорного распределения, нам
необходимо имитировать данную цепь Маркова согласно описанному выше алгоритму и
+N ·K
получить {h0j }M
элементов данной цепи. Затем необходимо отсеять первые M элеменl=j

тов и из оставшихся выбрать каждый N -й элемент. Полученное множество будет являться
0
аппроксимацией искомой выборки {h0j }N
j=1 [8], а так как любая история выполнений h од-

нозначно определяет значение всех случайных величин в вероятностной программе, то и
аппроксимацией выборки из P (T | X = x̂).
K выбирается больше единицы, чтобы исключить автокорреляцию h0i , которая естественным образом возникает в цепи Маркова. M выбирается достаточно больши́м, чтобы
независимо от выбора начальной точки h01 цепь успела «забыть» данный первоначальный
выбор; это особенно важно, когда начальное априорное распределение P (T, X) очень сильно отличается от апостериорного P (T | X = x̂). Какого-то общего правила для выбора
данных величин нет, и обычно они выбираются эмпирически.
1.2.5

Программная реализация статистического вывода

В предыдущем подпункте теоретически был описан алгоритм для обобщенного статистического вывода в вероятностных языках программирования.
Простая программная реализация впервые достаточно подробно была описана в
[16]:
1. Для инициализации h01 выбирается случайным образом путем выполнения вероятностной программы и фиксации наблюдаемых случайных величин в соответствии с имеющимися данными. При этом в памяти компьютера сохраняется
база данных активных случайных величин вместе с их значениями, а также сохраняется P (h01 ), которое подсчитывается во время первоначального выполнения
программы с учетом вероятности наблюдаемых случайных величин принять то
19

значение, которое они приняли. В [16] описывается схема адресации, которая
позволяет различать случайные величины между собой.
2. Затем для выбора последующего h0t+1 каждый раз из базы данных случайным
равномерным образом выбирается случайная величина r и она случайным образом варьируется в соответствии со своим локальным распределением предлагаемых переходов κ( · | r = r̂, h0t ), которое предварительно задается для всех
используемых примитивных (несоставных) случайных величин. Создается копия
базы данных, в которой значение случайной величины r заменяется на новое.
3. Вероятностная программа выполняется еще один раз, при этом случайные величины, для которых уже в базе данных имеется запись (согласно схеме адресации), принимают соответствующие старые значения, а r принимает свое новое
полученное значение.
4. В общем случае, если r влияет на поток выполнения программы, генерируются
значения для новых случайных величин, которые становятся активными. Для
этих случайных величин прежде в базе данных не было записей.
5. Также если r влияет на поток выполнения программы, то некоторые случайные
величины могут перестать быть активными. В данном случае соответствующие
им записи в базе данных будут невостребованы.
6. После выполнения программы второй раз можно считать, что мы получили h0∗ .
Также отметим, что при выполнении программы h0∗ мы получили и использовали
всё необходимые компоненты для подсчета P (h0∗ ) и Q(h0∗ | h0t ).
7. Вероятность обратного перехода Q(h0t | h0∗ ) может быть получена разными способами. Если в вероятностной программе есть только независимые случайные
величины (и нет перестановочных), то первые две компоненты Q(h0t | h0∗ ) могут
быть найдены тривиально, так как мы знаем количество случайных величин в
h0∗ и можем посчитать κ(r̂ | r = r∗ , h0∗ ). Для третьей компоненты нам нужно посчитать произведение вероятностей случайных величин, которые перестали быть
активными в h0∗ , то есть все невостребованные случайные величины в базе данных. В случае наличия и использования перестановочных случайных величин
мы можем имитировать ситуацию, что h0∗ является нашим старым состоянием,
а новое состояние мы получаем в точном соответствии с h0t [3].

20

8. Имея P (h0t ), P (h0∗ ), Q(h0∗ | h0t ) и Q(h0t | h0∗ ), мы подсчитываем коэффициент принятия α для алгоритма Метрополиса-Гастингса, и либо принимаем h0∗ с вероятностью α, что означает, что в следующий раз мы будем использовать уже новую
базу данных случайных величин со значением, соответствующими h0∗ , либо отклоняем с вероятностью (1 − α), что означает, что на следующем шаге h0t+1 мы
снова будем использовать старую базы данную от h0t .
9. Алгоритм повторяется с шага № 2.
1.3

Эффективность вывода

При использовании вероятностных языков программирования встает вопрос об эффективности статистического вывода, другими словами о том, как быстро мы генерируем
выборку желаемой точности из апостериорного распределения. Описанная в секции 1.2.5
программная реализация статистического вывода является неэффективной, так как происходит перевыполнение всей вероятностной программы, хотя вариация случайной величины r обычно имеет только локальный эффект.
Рассмотрим более подробно данную проблему на примере следующей вероятностной программы:
1 [ASSUME r a i n y −s e a s o n ( b e r n o u l l i 0 . 2 ) ]
2 [ASSUME c l o u d y
3

( b e r n o u l l i ( i f r a i n y −s e a s o n 0 . 8 0 . 3 ) ) ]

4 [ASSUME r a i n
5

( b e r n o u l l i ( i f cloudy 0.8 0 . 2 ) ) ]

6 [ASSUME s p r i n k l e r
7

( b e r n o u l l i ( i f cloudy 0.1 0 . 5 ) ) ]

8 [ASSUME wet−g r a s s
9
10
11

( bernoulli
( i f sprinkler ( i f rain 0.99 0.9)
( i f rain 0.9 0 . 0 1 ) ) ) ]
Данная вероятностная программа описывает статистически простую упрощенную
модель зависимости между сезоном, облачностью, дождем, работой разбрызгивателя и
состоянием травы (мокрая или нет) в какой-то день. Значение переменных следующее:

21

Рисунок 1 — Пример Байесовской сети. Из [20].
rainy-season: входит ли тот день в сезон дождей или нет?, cloudy: облачно в тот день
или нет?; rain: был ли дождь в тот день?; sprinkler: работал ли разбрызгиватель в тот
день или нет?; wet-grass: была ли трава мокрой в тот день? Данная модель может быть
представлена с помощью Байесовской сети доверия (см. рис. 1) и таблицами условных
вероятностей. Отметим, что любая Байесовская сеть доверия может быть представлена в
виде вероятностной программы на языке Church/Venture/Anglican, но не любая вероятностная программа может быть представлена Байесовской сетью.
В случае, если во время очередной итерации алгоритма Метрополиса-Гастингса в
качестве варьируемой случайной величины выбрана случайная величина rain (соответствующий узел на рис. 1 выделен самым темным цветом), для подсчета коэффициента
принятия достаточно лишь рассмотреть значения и вероятности при данных значениях
узлов «Дождь» и «Трава мокрая». Все остальные значения и их вероятности останутся
прежними. Иллюстрацию распространения возмущений в связи с вариацией случайной
величины «Дождь» см. на рис. 2.
В подобном простом примере этот факт не играет большой роли, так как отношение
случайных величин, требующих «переучета», к общему количеству активных случайных
величин невелико, однако при большом объеме данных это играет решающую роль при
выполнении статистического вывода во многих моделях при помощи программных реализаций вероятностных языков программирования. Например, в скрытых марковских
моделях [17, 2, 4] или в латентном размещении Дирихле [21, 20].
Например, в простой скрытой Марковской модели есть N скрытых и N наблюдаемых величин. При реализации статистического вывода алгоритмом, описанным в сек22

Рисунок 2 — Байесовская сеть в виде «отпечатка» выполнения вероятностной программы.
Из [20].
ции 1.2.5, одна итерация алгоритма Метрополиса-Гастингса имеет сложность (по времени)
O(N ), хотя желаемая и возможная сложность O(1) или по крайней мере O(log N ).
Для достижения желаемой сложности в каждой истории выполнений необходимо
отслеживать зависимости между случайными величинами. Следует отметить, что эти зависимости в общем случае могут изменяться в вероятностных программах. Например, в
следующей программе
1 [ASSUME a ( b e r n o u l l i 0 . 5 ) ]
2 [ASSUME b (gamma 1 1 ) ]
3 [ASSUME c ( normal ( i f a b 3 . 0 ) 1 ) ]
случайная величина c, т.о. (normal ...), зависит от случайной величины b, т.о. (gamma
...), не всегда, а только если случайная величина a принимает значение «ИСТИНА».
Описание структур данных и алгоритмов, необходимых для отслеживания зависимостей в режиме реального времени, были предварительно приведено в [20] и [22], а затем
более обширно и подробно в [3]. При использовании данных структур данных и алгоритмов временная сложность одной итерации алгоритма Метрополиса-Гастингса в простой
скрытой Марковской модели равна O(log N ), при этом если использовать упорядоченный

23

перебор случайных величин, то временная сложность снизиться до O(1), так как логарифмический фактор появляется в связи с необходимостью выбирать случайным образом
следующий узел (т.е. случайную величину) для вариации.
В простой скрытой Марковской модели количество «активных» случайных величин постоянно и вид зависимости тривиален (от каждой скрытой случайной величины
напрямую зависит только одна наблюдаемая случайная величина и только следующая
скрытая случайная величина). В более сложных моделях виды зависимостей более изощрены и количество случайных величин меняется от одной истории выполнений к другой. С
другой стороны, для большого количества порождающих моделей, используемых в настоящее время в машинном обучении, необходимо, чтобы время на одну итерацию оставалось
постоянным или росло хотя бы логарифмически вместе с линейным ростом количества
наблюдений, иначе статистический вывод будет невозможен за разумное время.
На рис. 3 показаны результаты применения алгоритмов и структур данных, описанных в [20, 22, 3]. При использовании старого подхода, описанного в [16] (см. секцию 1.2.5)
время на N итераций алгоритма Метрополиса-Гастингса росло квадратично с линейным
ростом размерности модели (N скрытых и N наблюдаемых случайных величин), а при
использовании предлагаемого подхода время растет квазилинейно.
На этом же рисунке видно, что благодаря локализации выполнения одной итерации Метрополиса-Гастингса стало возможным проводить приближенный статистический
вывод параллельно.
Следует отметить, что пространственная сложность алгоритма (объем памяти,
необходимый для его применения) равна, грубо говоря, O(N + K), где O(K) — стоимость хранения в памяти зависимостей между случайными величинами. В общем случае
эта величина может быть достаточно большой, и в работе [22] описаны возможности более
эффективного расходования памяти.
1.4

Порождающее вероятностное программирование в распознавании образов

Как было отмечено в начале, цель вероятностного программирования — облегчить задачу
моделирования порождающих вероятностных моделей и проведения вывода в них. Примером иллюстрации успешного предварительного применения вероятностного программирования может служить работа [5], в которой вероятностное программирование использует24

Среднее время на N итераций
алгоритма Метрополиса-Гастингса

160
140
120

Статистический вывод без оптимизации (красный)
Предложенные алгоритмы и структуры данных (зел.)
—//—, 2 потока (син.)
—//—, 4 потока (гол.)

100
80
60
40
20
0

50

100

150

200

250

Количество наблюдений N в скрытой модели Маркова

Рисунок 3 — Эффективность статистической вывода при его различных реализациях:
красный график соответствует простейшему алгоритму, описанному в 1.2.5, при котором
каждый раз происходит перевыполнение всей вероятностной программы. Зеленый график соответствует использованию новых алгоритмов и структур данных, что позволяет
производить статистический вывод асимптотически более эффективно. Синий и голубой
график показывают, что применение предлагаемого подхода также позволяет производить
статистический вывод параллельно, по крайней мере приближенно. Рис. из [20].
ся для моделирования вероятностной модели находящихся на изображении объектов и их
взаимодействия между собой. Байесовский подход к интерпретации изображений путем
задания априорного распределения на расположение объектов и на связи между ними
был предложен задолго до появления рассматриваемых языков программирования, однако именно с их появлением исследование и осуществление данного подхода стало проще,
так как вероятностные языки программирования позволяют композиционно и компактно
представлять вероятностные порождающие модели и проводить статистический вывод в
них.
В работе рассматриваются два примера: проблема графической CAPTCHA [23] —
«компьютерного теста, используемого для определения, кем является пользователь системы: человеком или компьютером», знакомого почти что каждому пользователю интернета; и проблема нахождения на изображении с камеры автомобиля дороги, разделительной
полосы, и левого и правого оврагов. Полученные результаты по своей эффективности на
рассматриваемых простых примерах не уступают другим современным подходам к решению этих задач, однако представление, моделирование и вывод проще осуществляется с

25

помощью вероятностного программирования.
1.5

О различиях между Church, Venture, Anglican

Целью данной главы, очевидно, не было подробное описание этих вероятностных языков, а
только краткое введение в них. Заинтересованному читателю мы можем порекомендовать
продолжить свое знакомство с вероятностным языком Church с [2, 24, 25, 16], Venture —
[3], Anglican — [26, 4, 27].
Хотя многое объединяет эти вероятностные языки, в некоторых принципиальных
вещах они различаются. Например, Venture на данный момент больше позиционируется
как универсальная платформа, включающая в себя разные виды алгоритмов и методов
статистического вывода с запроектированной возможностью добавлять новые с помощью
использования базовых компонент и методов. В рамках работы над Anglican развиваются методы обобщенного вывода с использованием методов фильтрации частиц. Church, с
другой стороны, позволяет производить статистический запрос внутри другого статистического запроса, что Venture и Anglican пока делать не могут.

26

2

Автоматическая генерация вероятностных программ

При использовании полных по Тьюрингу вероятностных языков программирования,
включающих в себя функции высших порядков3 , которыми в том числе являются языки
Church, Venture и Anglican, вероятностная программа одновременно является и порождающей моделью, и записанной процедурой для генерации элементов выборки из этой модели
путем выполнения исходного кода данной процедуры. Любая процедура в вероятностном
программировании является формально программным кодом, который описывает процесс
генерации элемента выборки при заданных аргументах данной функции. Таким образом,
процедуры вероятностных программ являются конструктивным способом описания условных распределений.
Полные по Тьюрингу и допускающие функции высших порядков вероятностные
языки программирования открывают возможность проведения вывода исходного текста
самих вероятностных программ, если задано априорное распределение на множестве исходного текста, с помощью операторов eval и apply. Грубо говоря, необходимо представить
вероятностную порождающую мета-модель, которая будет генерировать вероятностные
модели в виде исходного кода вероятностных программ.
Данная глава основана на работе [6], которая включает в себя первые предварительные результаты по этой амбициозной задаче вывода самих порождающих вероятностных
моделей при наличии какой-либо информации об искомом распределении, которое определяется искомой вероятностной программой. Отметим, что статистический вывод в пространстве исходного кода сложен и нет какого-то простого подхода как к построению
вероятностных порождающих моделей исходного кода [28], так и выводу в них [29].
В рамках нашей предварительной работы мы поставили задачу найти с помощью
статистического вывода и предлагаемой нами порождающей вероятностной мета-модели
такие вероятностные программы, которые будут генерировать выборку элементов, схожую по своим статистическим характеристикам с каким-то заранее заданным распределением.
Эта задача интересна сама по себе, так как нахождение эффективных алгоритмов
генерации (моделирования) случайных величин — нетривиальная задача и для людей3 Функциями (процедурами) высших порядков называются функции, аргументами или значениями которых могут быть
другие функции.

27

ученых, которой они занимаются на протяжении десятков лет [30, 31, 32].
Наши предварительные результаты показывают, что подобный автоматизированный вывод порождающих вероятностных моделей в виде вероятностных программ действительно возможен. В частности, мы приводим результаты успешного эксперимента, в
рамках которого мы автоматически нашли обобщенную вероятностную программу, которая подлинно (не приближенно) генерирует случайные величины, распределенные по
закону Бернулли с произвольным параметром p.
2.1

Обзор литературы

Рассматриваемые нами идеи относятся к разным областям, в том числе к автоматизации
процесса программирования [33, 34], индуктивному программированию [33, 35, 36, 37, 38,
39], автоматическому моделированию [40, 41], компьютерному определению и представлению плотности распределений. Подходы к решению проблемы включают статистический
вывод, поиск и оптимизацию, в том числе эволюционные алгоритмы и в особенности генетическое программирование.
Идеи и методы использования вероятностного программирования для изучения и
автоматизированного представления вероятностных моделей предлагались и ранее [2, 42,
43, 3].
Насколько автору известно, описываемый подход к порождению вероятностных
программам мета-вероятностной программой ранее не рассматривался в качестве отдельной проблемы достаточно подробно, хотя первые шаги в исследовании и формулировке
проблемы были сделаны в работах [42, 3].
2.2

Описание подхода

Наш подход может быть описан в рамках приближенных Байесовских вычислений [44]
с использованием метода Монте-Карло по схеме цепей Маркова с выбором в качестве
искомого апостериорного распределения
π(X |X̂ )p(X̂ |T )p(T ),

(1)

где π(X |X̂ ) — вероятностная мера, измеряющая расстояние между значения статистик,
вычисленных соответственно на выборке X из искомого распределения и выборке X̂ , по28

[ASSUME program−text ( p r o d u c t i o n s ‘ ( ) ‘ r e a l ) ]
[ASSUME program ( eval ( l i s t ‘ lambda ‘ ( ) program−text ) ) ]
[ASSUME samples ( apply−n−times program 10000 ’ ( ) ) ]
[OBSERVE ( normal ( mean samples ) n o i s e − l e v e l ) 0 . 0 ]
[OBSERVE ( normal ( v a r i a n c e samples ) n o i s e − l e v e l ) 1 . 0 ]
[OBSERVE ( normal ( sk ewn ess samples ) n o i s e − l e v e l ) 0 . 0 ]
[OBSERVE ( normal ( k u r t o s i s samples ) n o i s e − l e v e l ) 0 . 0 ]
[ PREDICT program−text ]
[ PREDICT ( apply−n−times . . . ( program ) ) ]
Рисунок 4 — Вероятностная программа для вывода исходного кода вероятностной программы для генерации случайных чисел, распределенных по закону стандартного нормального распределения Normal(0, 1).
лученной путем выполнения исходного кода вероятностной программы-кандидата.
Пусть существует распределение F с параметрами θ, вероятностную программу
для генерации элементов выборки из которого мы хотим вывести. Пусть X = {xi }Ii=1 , xi ∼
F (·|θ) будет выборкой из I элементов из распределения F при каком-то фиксированном
значении параметров θ. Рассмотрим задачу вывода исходного кода вероятностной программы T , которая при ее неоднократном выполнении J раз сгенерирует выборку элементов X̂ = {xˆj }Jj=1 , xˆj ∼ T (·), статистически близких к распределению F при заданных
параметрах θ.
Пусть s будет статистикой выборки, и тогда s(X ) и s(X̂ ) — значение этой статистики на элементах выборки из F и из T (·) соответственно. Пусть вероятностная мера
d(s(X ), s(X̂ )) = π(X |X̂ ), для простоты ненормированная, принимает бо́льшие значения,
когда s(X ) ≈ s(X̂ ). d можно интерпретировать как расстояние или штраф.
Мы используем вероятностное программирование для представления мета-модели,
порождающей другие вероятностные программы в виде исходного текста, и для проведения статистического вывода в пространстве искомых вероятностных моделей. Для статистического вывода мы использовали программную реализацию вероятностного языка программирования Anglican, которая поддерживает [4] статистический вывод методом частиц
и методом Метрополиса-Гастингса по методу Монте-Карло по схеме цепей Маркова.
Вероятностная мета-модель представлена на рис. 4, где на первой строке мы устанавливаем соответствие между T и переменной program-text, которая будет содержать
один сгенерированный элемент из распределения на исходный код P (T ), определенное
априорно через порождающую процедуру production с помощью адаптивной грамматики по типу [45] (см. подробнее в разделе 2.3).
29

Мы не указываем здесь θ, так как задача вывода в данном случае найти вероятностную программу, генерирующую элементы из стандартного нормального распределения.
Переменная samples на второй строке представляет описанную выше выборку X̂ из вероятностной программы-кандидата, и в этом примере J = 10000.
s и d вычисляются на следующих четырех строках вероятностной программы, где
статистика s определяется как четырехмерный вектор, включающий в себя соответственно
выборочные среднее, дисперсию, коэффициент асимметрии и коэффициент эксцесса выборки элементов из распределения, определенного вероятностной программой, полученной
из распределения T . Мера расстояния d определяется через плотность многомерного нормального распределения со средними [0.0, 1.0, 0.0, 0.0]T и диагональной ковариационной
матрицей σ 2 I. Отметим, что это означает, что мы ищем вероятностные программы, результат выполнения которых определяет распределение со средним равным 0, дисперсией
— 1, коэффициентами асимметрии и эксцесса равными 0, и мы «штрафуем» отклонения
от этих значений с помощью квадратичной экспоненциальной функции потерь с коэффициентом

1
,
σ2

где σ определена как noise-level. Эта функция потерь представляет собой

функцию плотности нормального распределения.
Данный пример служит хорошей иллюстрацией основных особенностей нашего подхода. Для поиска в виде вероятностной программы генератора случайных чисел из стандартного нормального распределения мы используем аналитическую информацию s(F ) о
стандартном нормальном распределении при вычислении расстояния между статистиками d(s(F ), s(X̂ )). Существует по крайней мере три различных ситуации, включая данную,
в которых s и d могут вычисляться различными способами:
1. В рамках первой ситуации мы ищем вероятностную программу, которая эффективно генерирует элементы из аналитически известных распределений. Под эффективностью в данном случае можно понимать вычислительную временную
сложность и среднее количество использованной энтропии для генерации элемента выборки в среднем. Практически всегда в этой ситуации и при подобной
постановке задачи статистики распределения F известны аналитически.
2. Вторая ситуация возникает, когда мы можем генерировать только элементы выборки из F . Например, подобная ситуация возникает, когда мы генерируем элементы из апостериорного распределения с помощью «дорогостоящего» (вычис-

30

[ASSUME program−text ( p r o d u c t i o n s ‘ ( r e a l ) ‘ b o o l ) ]
[ASSUME program ( eval ( l i s t ‘ lambda ‘ ( ) program−text ) ) ]
[ASSUME J 1 0 0 ]
[ASSUME samples−1 ( apply−n−times program J ’ ( 0 . 5 ) ) ]
[OBSERVE ( f l i p ( G−test−p−value
samples−1 ‘ B e r n o u l l i ( l i s t 0 . 5 ) ) ) t r u e ]
...
[ASSUME samples−N ( apply−n−times program J ’ ( 0 . 7 ) ) ]
[OBSERVE ( f l i p ( G−test−p−value
samples−N ‘ B e r n o u l l i ( l i s t 0 . 7 ) ) ) t r u e ]
[ PREDICT program−text ]
[ PREDICT ( apply−n−times program J ’ ( 0 . 3 ) ) ]
Рисунок 5 — Вероятностная программа для вывода исходного кода вероятностной
программы, генерирующей случайные числа, распределенные по закону Бернулли
Bernoulli(θ). На предпоследней строчке выводится текст вероятностной программыкандидата. На последней строчке вероятностная программа-кандидат выполняется J раз
для генерации выборки из J элементов при параметризации p = 0.3, причем предыдущие
θn (т.о. тренировочные значения параметров) не содержали p = 0.3.
лительно) метода Монте-Карло по схеме Марковских цепей, и мы заинтересованы получить другое представление апостериорного распределения в виде вероятностной программы, чье априорное распределение будет точно или хотя бы
приблизительно совпадать с искомым апостериорным.
3. При третьей ситуации нам заранее дана фиксированного размера выборка из
F , и мы хотим найти вероятностную программу, которая позволит генерировать
выборку произвольного размера из F в дальнейшем. Мы начали этот раздел с
постановки задачи именно в рамках третьей ситуации.
Следует отметить, что возможная польза представления генератора случайных чисел потенциально не только в том, чтобы эффективно и быстро генерировать выборку
из желаемого распределения, но и чтобы получить формальное представление желаемого
распределения или его приближения в виде исходного кода вероятностной программы, то
есть в виде формальной сущности, дальнейшие действия с которой можно производить; в
том числе проводить анализ над выведенным исходным кодом и использовать найденные
блоки исходного кода для решения других задач.
Рис. 5 иллюстрирует другое важное обобщение решение задачи, сформулированной
в самом начале с апостериорным распределением (1).
При постановке задачи на вывод генератора случайных чисел, распределенных согласно стандартному нормальному распределению, мы не параметризовали искомое рас-

31

[ASSUME box−muller−normal
( lambda ( mean s t d )
(+ mean ( ∗ s t d
( ∗ ( cos ( ∗ 2 ( ∗ 3 . 1 4 1 5 9
( uniform−continuous 0 . 0 1 . 0 ) ) ) )
( sqrt ( ∗ −2
( log ( u n i f o r m − c o n t i n u o u s 0 . 0 1 . 0 )
)))))))]

[ASSUME p o i s s o n ( lambda ( r a t e )
( b e g i n ( d e f i n e L ( exp ( ∗ −1 r a t e ) ) )
( d e f i n e i n n e r − l o o p ( lambda ( k p )
( i f (< p L) ( dec k )
( begin ( d e f i n e u
( uniform−continuous 0 1 ) )
( inner−loop ( inc k ) (∗ p u ) ) ) ) ) )
( inner−loop 1 ( uniform−continuous 0 1 ) ) ) ) ]

Рисунок 6 — Найденные и записанные людьми исходные коды вероятностных программ
для (слева) общего нормального распределения Normal(µ, σ) [31] и (справа) распределения
Пуассона Poisson(λ) [46]. Эти исходные коды входят в собранный нами корпус, с помощью
которого мы определяем априорные вероятности для наших порождающих правил путем
подсчета количества встречающихся констант, процедур разных видов и т.д.
пределение никаким образом, так как у стандартного нормального распределения нет
параметров, т.е. θ = ∅. В общем случае распределения, которые мы хотим представить
в виде вероятностных программ, имеют нетривиальную параметризацию, и представляют собой по сути семейство распределений. Мы хотим найти вероятностную программу,
входные аргументы которой как раз бы и являлись параметрами искомого распределения,
таким образом, эта вероятностная программа позволяла бы генерировать случайные величины из всего семейства. Для наглядности рассмотрим алгоритм генерации случайных
чисел, распределенных по закону нормального распределения, с помощью преобразования
Бокса-Мюллера, представленный в виде вероятностной программы на рис. 6.
Эта вероятностная процедура параметризована двумя параметрами: средним и
стандартным отклонением. Для постановки обобщенной задачи вывода параметризованных вероятностных программам нам необходимо изменить наше искомое апостериорное
распределение, включив в него параметр θ, параметризующий искомое распределение
F:
π(X |X̂ , θ)p(X̂ |T , θ)p(T |θ)p(θ).
Нам бы хотелось вывести вероятностную программу, которая смогла бы обобщить
все возможные значения параметра θ. С допущением, что если мы выберем конечное число N различных параметризаций θ̂i , мы получим обобщение всего семейства распределений в виде вероятностной программы, мы формулируем нашу задачу в виде следующего
приближения с использованием приближенных Байесовских вычислений в рамках метода

32

Монте-Карло по схеме цепей Маркова:
Z
N
1 X
π(Xn |X̂n , θn )p(X̂n |T , θn )p(T |θn ) ≈ π(X |X̂ , θ)p(X̂ |T , θ)p(T |θ)p(θ)dθ.
N n=1
Вероятностная программа для поиска параметризованной вероятностной программы, генерирующей случайные числа, распределенных по закону Бернулли Bernoulli(θ),
представлена на рис. 5, и наглядным образом иллюстрирует применение данного допущения и приближения. При выбранных N различных параметризациях параметра p распределения Бернулли мы каждый раз генерируем J элементов из вероятностной программыкандидата, аккумулируя расстояние (штраф) между искомым распределением и полученным распределением, представляющим вероятностную программу-кандидата. В каждом
конкретном случае для θi мы высчитываем: 1) расстояние с использованием статистики G-теста (более «современный» аналог критерия согласия Пирсона и соответствующей
статистики) в виде

Gn = 2

X
i∈0,1

#[X̂n = i]

#[X̂n = i]ln

θni (1 − θn )(1−i) · |X̂n |

!
,

где #[X̂n = i] — количество элементов выборки из X̂n , принимающих значение i; 2) а
также соответствующее p-значение с нулевой гипотезой, утверждающей, что элементы
выборки X̂n являются и элементами выборки из распределения Bernoulli(θn ). Так как
распределение статистики G-теста приблизительно распределено по закону распределения
хи-квадрат, т.е. G ∼ χ2 (1) в нашем примере, мы можем представить и находить расстояние
d в данном случае путем вычисления вероятности ложного отклонения нулевой гипотезы
H0 : X̂n ∼ Bernoulli(θn ). Ложное отклонение нулевой гипотезы эквивалентно успеху в
проведении испытания Бернулли с вероятностью успеха равной p-значению.
2.3

Грамматика и порождающие правила

С учетом наличия в нашем распоряжении выразительного вероятностного языка программирования, допускающего функции высших порядков и полного по Тьюрингу, наше
априорное распределение об исходном коде искомых вероятностных программ также достаточно выразительно. В общих чертах оно схоже с адаптивными грамматиками [45],
используемыми в [29], но имеет отличия, в частности связанные с созданием сред с ло-

33

кальными переменными. В виде псевдокода наше априорное распределение может быть
представлено следующим образом (символ → означает «может перейти в»):
1. Выражение exprtype |env → в имя переменной, случайно выбираемой из среды
переменных env с типом type.
2. Выражение exprtype |env → в случайную константу типа type. Константы различных типов (целочисленные, вещественные и т.д.) генерируются из отдельного для
каждого типа type процесса Дирихле4 DPtype (Htype , α), где базовое распределение Htype само по себе в общем случае являются смесью нескольких распределений. Например, для констант вещественного типа мы используем смесь нормального распределения (normal 0 10), равномерного непрерывного распределения
(uniform-continuous -100 100) и равномерного дискретного распределения из
множества {−1, 0, 1, . . .}.
3. Выражение exprtype |env → (proceduretype exprarg_1_type ... exprarg_N _type ), где процедура procedure является случайно выбираемой примитивной детерминированной или стохастической (не составной) процедурой, определенной заранее в глобальной среде, с типом возвращаемого значения type.
4. Выражение exprtype |env → (compound_proceduretype exprarg_1_type ... exprarg_N _type ),
где compound_proceduretype является составной процедурой, также генерируемой
в соответствии с процессом Дирихле DPtype (Gtype , β), где базовое распределение
Gtype случайным образом генерирует составную процедуру с типом возвращаемого значения type и количеством входных аргументов, распределенным
по закону Пуассона, где каждый входной аргумент имеет свой произвольный
тип. Тело body составной процедуры генерируется также случайным образом
согласно этим же порождающим правилам и грамматике, но с учетом введение
локальной среды, включающей в себя входные аргументы процедуры.
5. Выражение exprtype |env → (let [(gensym) exprreal ] exprtype |env ∪ gensym)), где
env ∪ gensym означает среду, дополненную переменной с именем (gensym) и со
значением, вычисляемым в соответствии с генерируемым случайным образом
выражением согласно этих же порождающих правил.
6. Выражение exprtype |env → (if (exprbool ) exprtype exprtype ).
7. Выражение exprtype |env → (recur exprarg_1_type ... exprarg_M _type ), таким образом
4 Не

нужно путать с распределением Дирихле.

34

рекурсивный вызов текущей составной процедуры.
Во избежание вычислительных ошибок во время выполнения сгенерированных процедур мы заменяем примитивные функции их «защищенными» аналогами, например
log(a) заменяется на safe-log(a), причем последний возвращает 0 если a < 0; или например uniform-continuous заменяется safe-uc(a, b), которая в случае если a > b меняет
аргументы местами, а также возвращает просто a если аргументы равны: a = b.
Множество типов, которые мы использовали в рамках наших экспериментов, включало вещественные и булевы типы, а общее множество примитивных процедур, включенных нами в глобальную среду, включало в себя такие функции как +, −, *, safe-div,
safe-uc, safe-normal.
2.4

Вероятности использования порождающих правил

Для задания априорных вероятностей порождающих правил, то есть вероятностей, с которой каждое из правил будет применяться в случае возможности его применения, мы
вручную составили небольшой корпус вероятностных программ, которые повторяют найденные учеными [32] алгоритмы генераторов случайных чисел. Примеры таких программ
представлены на рис. 6. Заметим, что все они требуют наличия только одной стохастической процедуры, а именно uniform-continuous, так что мы включили только ее в глобальную среду с положительной вероятностью для экспериментов, описанных в 2.5.2.
Используя данный корпус, мы вычислили априорные вероятности каждого порождающего правила, при этом при выводе вероятностной программы для генерации случайных величин из искомого распределения F (например, распределения Бернулли), мы
исключали из корпуса все элементы, которые генерируют случайные величины согласно
закону распределения F . После этого вероятности использования порождающих правил
были «смягчены» с помощью распределения Дирихле. В будущем можно использовать более обширные корпусы вероятностных программ, примером зарождающегося подобного
корпуса может служить [25].
2.5

Эксперименты

Эксперименты были разработаны таким образом, чтобы проиллюстрировать все три вида
возможных ситуаций, описанных ранее. Но в начале мы проиллюстрируем выразительность нашего априорного распределения исходного кода вероятностных программ. После
35

−4
Y

−3

−2

0.6
0.5
P̂ (Y )

P̂ (Y )
−5

0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.00

−2.5

−2.0

−1.5
Y

−1.0

−0.5

0.0

0.0

1.0
0.8
0.6
0.4
0.2
−0.5

0.0
Y

0.5

1.0

0.0

0.3
0.1

−3.0

1.2

−1.0

0.4
0.2

P̂ (Y )

1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

−6

P̂ (Y )

P̂ (Y )
P̂ (Y )

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

−14.56636

Y

−7.28318

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

−4

−2

0
Y

2

4

−2.0

−1.5

−1.0
Y

−0.5

0.0

Рисунок 7 — Гистограммы выборок из некоторых порожденных вероятностных программ
из их априорного распределения. По форме распределений видно, что наши порождающие правила достаточно обширны, чтобы соответствующие порождающим вероятностным
программам распределения были нетривиальны.
этого мы опишем постановку и результаты экспериментов в рамках использования нашего
подхода во всех трех различных ситуациях возможности вычисления расстояния d.
2.5.1

Выборки из сгенерированных вероятностных программ

Для иллюстрации выразительности нашего априорного распределения исходных кодов
вероятностных программ мы приводим выборки из случайным образом сгенерированных
вероятностных программ из их априорного распределения. Эти шесть выборок в виде гистограмм из шести различных автоматически сгенерированных вероятностных программам расположены на рис. 7.
Из рисунка видно, что разные случайным образом сгенерированные вероятностные
программы определяют в общем случае достаточно различные структурно распределения.
В частности, можно заметить разнообразие в носителе, дисперсии и количестве «режимов»
(т.е. островков носителя с высокой вероятностью по отношению к очень низковероятностным пространствам между островками).
2.5.2

Вывод вероятностных программ, определяющих распределения, совпадающие или
приближающие классические одномерные распределения

Как уже было отмечено, для всех классических одномерных распределений существуют
алгоритмы, позволяющие и точно, и приближенно генерировать любое количество элементов из данных параметризованных распределений, и они были аналитически выведены
учеными. Эти алгоритмы в том числе, очевидно, могут быть записаны как вероятностные
программы в виде их исходного кода.
36

Мы провели серию экспериментов, чтобы проверить возможность автоматического вывода вероятностных программ, генерирующих выборки из классических одномерных
распределений, а именно из распределения Бернулли Bernoulli(p), распределения Пуассона
Poisson(λ), Гамма-распределения Gamma(a, 1.0), Бета-распределения Beta(a, 1), стандартного нормального распределения Normal(0, 1), и «общего» нормального распределения
Normal(µ, σ).
Для каждого семейства распределений мы проводили статистический вывод,
выбирая в качестве целевого апостериорного распределения частное распределение
π(X |X̂ )p(X̂ |T )p(T ), маргинализированное по параметру θ. Для приближения в каждом
случае мы выбирали малое множество значений параметров {θ1 , . . . , θN } и определяли апостериорное распределение ограничениями по p−значению (для распределения Бернулли)
и по близости моментам к ожидаемым (для всех других рассматриваемых распределений).
Еще раз отметим, что при задании априорного распределения на порождающие правила
в каждом конкретном случае из корпуса исключались все вероятностные программы, относящиеся к искомому распределению.
Образцы гистограмм выборок из лучших найденных в результате вывода вероятностных программ представлены на рис. 8, где под лучшими мы понимаем вероятностные
программы с наименьшим расстоянием d на тренировочных значениях параметров и соответствующих значениях моментов или максимального p-значения.
Стоит особенно отметить результат эксперимента с распределением Бернулли, в
рамках которого был найден исходный код вероятностной программы, статистически подлинно (точно) генерирующей выборку из всего семейства распределения Бернулли с параметром p. Найденная вероятностная программа представлена на рис. 9.
На рис. 10 представлен выведенный исходный текст вероятностной программы для
генерации элементов из Гамма-распределения Gamma(a, 1), параметризованного параметром a.
2.5.3

Обобщение произвольных эмпирических данных с помощью порождающего вероятностного программирования

Мы также проверили наш метод на выводе порождающих моделей в виде вероятностных
программ для объяснения настоящих, не синтетических, данных, аналитическое распределение которых неизвестно.
37

1.2

0.05
0.04

0.6

0.03

0.4

0.02

0.2

0.01

0.0

0

0.00

1

P̂ (Y )

0.8

P̂ (Y )

P̂ (Y )

1.0

40

50

60

70

0

2

4

6
Y

90

100

110

8

10

12

14

0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.00

-2 -1 0 1 2 3 4 5 6 7 8 9 10 11
Y

2.0
1.5
P̂ (Y )

0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.00

80
Y

P̂ (Y )

P̂ (Y )

Y

0.45
0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
0.00

1.0
0.5

−4

−2

0
Y

2

4

0.0

0.0

0.5

1.0

1.5

Y

Рисунок 8 — (Зеленые сплошные линии) Гистограммы выборок из распределений, соответствующих вероятностным программам, имеющим «высокую» вероятность в апостериорном распределении при поиске порождающих моделей, соответствующих распределениям (слева направо, сверху вниз): Bernoulli(p), Normal(µ, σ), Poisson(λ), Gamma(a, 1.0),
Normal(0, 1), Beta(a, 1.0). (Синие пунктирные линии) Гистограммы выборок из «настоящих» распределений. Параметризация семейства распределений в каждом случае производилась другими значениями параметров, которые не входили в обучающее множество
значений {θn }. Отметим, что для семейства распределений Бернулли Bernoulli(p) была выведена вероятностная программа, статистически точно генерирующая элементы выборки
при любой параметризации (см. рис. 9). С другой стороны, не все распределения получили хорошее приближение найденными вероятностными программами за то время, что мы
проводили вывод, как, например, в случае Бета-распределение Beta(a, 1.0).

( lambda ( par s t a c k − i d ) ( i f (< ( uniform−continuous 0 . 0 1 . 0 ) par ) 1 . 0 0 . 0 ) )
( lambda ( par
( i f (< 1 . 0
( lambda ( par
( i f (< 1 . 0

stack−id )
( s a f e − s q r t ( s a f e − d i v par ( safe−uc par ( dec par ) ) ) ) ) 1 . 0 0 . 0 ) )
stack−id )
( safe−uc ( s a f e − s q r t par ) (+ par ( cos par ) ) ) ) 1 . 0 0 . 0 ) )

Рисунок 9 — (сверху) Написанный человеком исходный код генератора случайных чисел,
распределенных по закону Бернулли Bernoulli(p). (внизу, две программы) Выведенные исходные коды. Первая из двух выведенных вероятностных программ определяет настоящее
семейство распределений Бернулли Bernoulli(p), параметризованное p. Вторая программа
генерирует распределение, приближенное к распределению Бернулли, параметризованное
p.

38

( lambda ( par s t a c k − i d ) ( ∗ ( b e g i n ( d e f i n e sym0 0 . 0 )
( exp ( s a f e − u c −1.0 ( s a f e − s q r t ( s a f e − u c
( safe−div ( safe−uc 0 . 0 ( safe−uc 0 . 0 3 . 1 4 1 5 9 ) )
par ) (+ 1 . 0 ( s a f e − u c ( b e g i n ( d e f i n e sym2
( lambda ( v a r 1 v a r 2 s t a c k − i d ) ( dec v a r 2 ) ) )
( sym2 ( s a f e − u c −2.0 ( ∗ ( s a f e − u c 0 . 0 ( b e g i n
( d e f i n e sym4 ( s a f e − u c sym0 ( ∗ (+ ( b e g i n
( d e f i n e sym5 ( lambda ( v a r 1 v a r 2 s t a c k − i d )
( s a f e − d i v (+ ( s a f e − l o g ( dec 0 . 0 ) ) −1.0) v a r 1 ) ) )
( sym5 ( exp par ) 1 . 0 0 ) ) 1 . 0 ) 1 . 0 ) ) )
( i f (< ( s a f e − u c par sym4 ) 1 . 0 ) sym0
( s a f e − u c 0 . 0 − 1 . 0 ) ) ) ) sym0 ) )
( s a f e − d i v sym0 ( exp 1 . 0 ) ) 0 ) ) 0 . 0 ) ) ) ) ) ) ) par ) )

( lambda ( s t a c k − i d )
(∗ 2.0 (∗ (∗
( ∗ −1.0 ( s a f e − u c 0 . 0 2 . 0 ) )
( safe−uc ( safe−uc 4 . 0
(+ ( s a f e − l o g 2 . 0 ) − 1 . 0 ) )
(∗ ( safe−div 2.0
−55.61617747203855)
( i f (< ( s a f e − u c
( safe−uc
27.396810474207317
( s a f e − u c −1.0 2 . 0 ) )
2.0) 2.0)
4.0 −1.0)))) −1.0)))

0.20

0.25

0.04

0.15

0.20

0.03
0.02

0.10
0.05

0.01
0.00

P̂ (Y )

0.05

P̂ (Y )

P̂ (Y )

Рисунок 10 — Исходный код выведенных вероятностных программ для (слева) Гаммараспределения Gamma(a, 1) и (справа) для третьей эмпирической выборки индикаторов,
используемых для рассмотрения заявок на выдачу кредита. Выборки, сгенерированные с
помощью выполнения этих программ, расположены соответственно на рис. 8 и на рис. 11
(последняя гистограмма из трех). Для экономии места исходный код был сокращен, где
возможно; например, путем замены (* 1.0 0.0) на 0.0.

10

20

30

40

50

60

70

80

90

0.00

0.15
0.10
0.05

0

5

10

Y

15
Y

20

25

30

0.00
−5

0

5

10

15
Y

20

25

30

Рисунок 11 — Гистограммы (зеленые сплошные линии) выборок из распределений, определяемых найденными в процессе статистического вывода вероятностными программами
для аппроксимации эмпирических данных (синие пунктирные линии) трех вещественных
показателей из базы данных, используемой для анализа заявок на выдачу кредита.
Для данного эксперимента мы выбрали три набора данных признаков обращающихся в банк клиентов для получения кредита [47, 48], и производили вывод вероятностной
программы, используя для сравнения p-значения двухвыборочного теста КолмогороваСмирнова аналогично тому, как мы использовали G-тест для дискретных распределений. Гистограммы выборок из лучших найденных вероятностных программ в сравнении
с гистограммами истинных эмпирическими данными признаков клиентов приведены на
рис. 11. Пример выведенной вероятностной программы показан на рис. 10 (справа).
2.5.4

«Компиляция» вероятностных программ

Генерация выборок из апостериорного распределения методом Монте-Карло по схеме Марковских цепей, особенно в случае Байесовского вывода, обычно достаточно дорогостояща.
Под «компиляцией» вероятностных программ мы имеем в виду поиск вероятностных программ, априорное распределение которых согласовывалось бы точно или приблизительно
с искомым апостериорным распределением.
В вероятностном программировании, где в общем случае генерация выборки из же39

[ASSUME t h e t a ( b e t a 1 . 0 1 . 0 ) ]
[OBSERVE ( f l i p t h e t a ) True ]
[OBSERVE ( f l i p t h e t a ) True ]
[OBSERVE ( f l i p t h e t a ) True ]
[OBSERVE ( f l i p t h e t a ) True ]
[ PREDICT t h e t a ]
[ PREDICT ( f l i p t h e t a ) ]

[ASSUME t h e t a ( s a f e − b e t a 4 . 4 4 0 1 . 0 ) ]
[ASSUME t h e t a ( s a f e − s q r t
( safe−beta ( safe−log 11.602) 1 . 0 ) ) ]
[ASSUME t h e t a ( s a f e − b e t a ( s a f e − s q r t 2 7 . 8 1 0 ) 1 . 0 ) ]
[ASSUME t h e t a ( b e t a 5 . 0 1 . 0 ) ]
[ PREDICT t h e t a ]

Рисунок 12 — (слева) Вероятностная модель бета-биномиального распределения (в несжатой форме) в виде вероятностной программы. В рамках эксперимента была поставлена
задача найти формализацию апостериорного распределения скрытого параметра θ в виде другой вероятностной программы, чье априорное распределение будет совпадать или
приближать данное апостериорное. (справа, сверху) Найденные вероятностные программы (т.е., результат «компиляции» вероятностных программ), априорное распределение
которых приближает заданное апостериорное распределение. (справа, внизу) Записанный
человеком исходный код вероятностной программы, чье априорное распределение совпадает с априорным распределением. В данном конкретном случае данная вероятностная
программа может быть просто выведена аналитически, так как Бета-распределение является сопряженным к биномиальному, но в общем случае это нетривиальная задача.
лаемого апостериорного распределения путем использования методом Монте-Карло или
аналогичным им является единственным доступным средством за разумное время, эта
проблема стоит еще острее.
В качестве самых предварительных результатов мы провели эксперимент, в рамках которого наша априорная модель была моделью бета-биномиального распределения (в
несжатой форме) с априорным распределением на скрытый параметр θ ∼ Beta(1.0, 1.0),
и использовали алгоритм Метрополиса-Гастингса для получения выборки X̂ из апостериорного распределения скрытого параметра θ с учетом проведения четырех успешных
испытаний по схеме Бернулли Bernoulli(θ). Соответствующая вероятностная программа
представлена на рис. 12.
Затем мы использовали наш подход к выводу вероятностных программ, априорное
распределение которых будет статистически схожим с полученной выборкой из желаемого апостериорного распределения. Примеры найденных вероятностных программ даны
на рис. 12. В данном конкретном случае мы можем аналитически найти и записать в
виде вероятностной программы апостериорное распределение, равное Beta(5.0, 1.0). Таким образом, полученные результаты показывают, что мы нашли хорошее приближение
апостериорного распределения.

40

2.6

Обобщение порождающего вероятностного программирования

Как было отмечено ранее, порождающие вероятностные модели определяют совместное
распределение (T, X), часто задаваемое сначала с помощью распределения (T ), а затем
с помощью условного распределения (X | T ). Порождающие вероятностные модели могут использоваться либо для генерации выборок напрямую из (T, X), либо в качестве
промежуточного этапа для нахождения условного распределения скрытых параметров
(T | X = x̂), обычно с помощью приближенных методов статистического вывода или оптимизации.
В рамках наших экспериментов, описанных выше, мы решали задачу вывода вероятностной программы, которая определяет вероятностную модель (в общем случае параметризованную), априорное распределение которой соответствует или приближает искомое распределение, заданное аналитически, или в виде выборки, или в виде дорогостоящего генератора элементов выборки.
На языке вероятностного программирования Church, поддерживающего проведение
статистического запроса внутри другого статистического запроса, данная задача может
быть сформулирована следующим образом (подобно тому, как мы ее формулировали в
Anglican/Venture, см. 4 и 5):
1 ( query
2

( d e f i n e program−text ( p r o d u c t i o n s INPUT−TYPES OUTPUT−TYPE) )

3

( d e f i n e program

4

( eval ( l i s t ‘ lambda ‘ ( arg1 . . . ar g N) program−text ) ) )

5

program−text

6

(and

7
8

( noisy−distance−equal
( get−statistics

9

( multiple−query

10

; No d e f i n e s .

11

( program arg1_1 . . . arg1_N )

12

true ))

13
14

expected−statistics1 )
...
41

( noisy−distance−equal

15

( get−statistics

16
17

( multiple−query

18

; No d e f i n e s .

19

( program argM_1 . . . argM_N)

20

true ))
expected−statisticsM )))

21

где (query defines expression predicate) — задание запроса в Church (см. [2]),
INPUT-TYPES — перечисление типов входных аргументов, OUTPUT-TYPE — тип выходного значения процедуры, noisy-distance-equal — функция сравнения статистик,
get-statistics — функция извлечения статистик из распределения, multiple-query —
аналог query, возвращающий несколько элементов из распределения вместо одного элемента, expected-statistics... — M значений статистик искомого распределения.
Автор считает, что более общо, для нахождения вероятностной модели, служащей
вспомогательным инструментом в нахождении условного апостериорного распределения
P (T | X = x̂), мы можем быть заинтересованы в постановке задачи следующим образом:
1 ( query
2

( d e f i n e l a t e n t − v a r i a b l e s − c r e a t o r ( l a t e n t − p r o d u c t i o n s LATENT−TYPES) )

3

( d e f i n e program−text

4
5
6

( p r o d u c t i o n s LATENT−TYPES INPUT−TYPES OUTPUT−TYPE) )
( d e f i n e program
( eval ( l i s t ‘ lambda ‘ ( arg1 . . . argN ) program−text ) ) )

7

( l i s t l a t e n t − v a r i a b l e s − c r e a t o r program−text )

8

(and

9

( noisy−latents−equal

10

( multiple−query

11

( define latent−variables ( latent−variables−creator ))

12

latent−variables

13

( noisy−outputs−equal

14

( program l a t e n t − v a r i a b l e s arg1_1 . . . arg1_N ) output1 ) )
42

15

expected−latents1 )

16

...

17

( noisy−latents−equal

18

( multiple−query

19

( define latent−variables ( latent−variables−creator ))

20

latent−variables

21

( noisy−outputs−equal
( program l a t e n t − v a r i a b l e s argM_1 . . . argM_N) outputM ) )

22
23

expected−latentsM ) ) )
где latent-productions — правила порождения объекта, представляющего скрытые параметры типов LATENT-TYPES, «лениво» создаваемые с помощью процедуры mem (см. [2]);
latent-variables — объект, представляющий скрытые параметры; arg..._... — M значений аргументов, output... — M значений наблюдаемых данных; expected-latents...
— M значений ожидаемых скрытых параметров.
Таким образом, в рамках данной формулировки мы ищем вероятностную программу, содержащую: скрытые переменные, определяемые latent-variables; наблюдаемые
переменные, определяемые выполнением program; и вывод производится с учетом M тренировочных троек
{ai , t̂i , x̂i },
где ai — дополнительная параметризация, x̂i — значения наблюдаемых переменных, t̂i —
ожидаемые значения скрытых переменных.
Следует отметить, что так как в настоящее время обобщенный статистический вывод с использованием вероятностных языков программирования производится достаточно
медленно, то решение практических задач таким образом пока не имеет смысла. С другой
стороны, с дальнейшим развитием [4, 3] обобщенных методов статистического вывода и
других решений (например, оптимизации или поиска), с использовании более продуманных алгоритмов программной реализации и структур данных [3, 27], со специализированным аппаратным обеспечением [49, 50], автоматизированная генерация и использование порождающих вероятностных моделей может успешно и продуктивно применяться на
практике.

43

ЗАКЛЮЧЕНИЕ
В первой части данной работы очень кратко и реферативно представлено введение
в вероятностное программирование на примере языков Church/Venture/Anglican, первое
подобное введение на русском языке, насколько известно автору.
Во второй части данной работы представлен подход к порождению вероятностных программ, обобщающих распределения, представленные либо в виде выборки, либо в
виде аналитического представления (например, в виде значения статистик), причем метавероятностная модель вероятностных моделей также записывается в виде вероятностной
программы. Используя данный подход, были получены предварительные положительные
результаты статистического вывода из апостериорного распределения искомых вероятностных программ методами Монте-Карло по схеме цепей Маркова, в том числе была автоматически выведена вероятностная программа, генерирующая элементы выборки всего
семейства распределений Бернулли.

44

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ
[1] Goodman Noah D. The principles and practice of probabilistic programming // ACM
SIGPLAN Notices / ACM. — Vol. 48. — 2013. — P. 399–402.
[2] Church: A language for generative models / Noah D. Goodman, Vikash K. Mansinghka,
Daniel M. Roy et al. — 2008. — P. 220–229.
[3] Mansinghka Vikash, Selsam Daniel, Perov Yura. Venture: a higher-order probabilistic
programming platform with programmable inference // arXiv preprint arXiv:1404.0099. —
2014.
[4] Wood Frank, van de Meent Jan Willem, Mansinghka Vikash. A new approach to
probabilistic programming inference // Proceedings of the 17th International conference
on Artificial Intelligence and Statistics. — 2014.
[5] Approximate Bayesian image interpretation using generative probabilistic graphics
programs / Vikash Mansinghka, Tejas D Kulkarni, Yura N Perov, Josh Tenenbaum //
Advances in Neural Information Processing Systems. — 2013. — P. 1520–1528.
[6] Perov Yura N, Wood Frank D. Learning probabilistic programs (статья отправлена на
конференцию). — 2014.
[7] Лопатников Леонид Исидорович. Экономико-математический словарь. — 5-е изд., перераб. и доп. — М.: Дело, 2003.
[8] Ветров Дмитрий Петрович, Кропотов Дмитрий Александрович. Байесовские методы машинного обучения (курс лекций). Методы Монте-Карло по схеме Марковских
цепей. — 2011. — URL: http://www.machinelearning.ru/wiki/images/6/6b/BMMO11_
10.pdf.
[9] Рассел Стюарт, Норвиг Питер. Искусственный интеллект: современный подход (2-е
издание). — «Вильямс», 2007.
[10] Bishop Christopher. Pattern Recognition and Machine Learning. — Springer, New York,
2006.
[11] Murphy Kevin P. Machine learning: a probabilistic perspective. — MIT Press, 2012.
[12] Ветров Дмитрий Петрович, Кропотов Дмитрий Александрович. Байесовские методы

45

машинного обучения (курс лекций). — 2013. — URL: http://www.machinelearning.
ru/.
[13] Сайт

о

вероятностном

программировании. —

2014. —

URL:

http://

probabilistic-programming.org/.
[14] Абельсон Харольд, Сассман Джеральд Джей. Структура и интерпретация компьютерных программ. — Добросвет, КДУ, 2010.
[15] Freer Cameron, Mansinghka Vikash, Roy Daniel. When are probabilistic programs
probably computationally tractable? // NIPS 2010 Workshop on Monte Carlo Methods
in Modern Applications. — 2010.
[16] Wingate David, Stuhlmueller Andreas, Goodman Noah D. Lightweight implementations of
probabilistic programming languages via transformational compilation. — 2011. — P. 131.
[17] Николенко Сергей Игоревич. Вероятностное обучение (курс лекций). — 2007. — URL:
http://logic.pdmi.ras.ru/~sergey/index.php?page=mlbayes.
[18] Chib Siddhartha, Greenberg Edward. Understanding the metropolis-hastings algorithm //
The American Statistician. — 1995. — Vol. 49, no. 4. — P. 327–335.
[19] Зубков А. М., Шуваев Д. В. Вычисление моментов комбинаторных статистик от перестановочных случайных величин // Дискретная математика. — 2005.
[20] Perov Yura, Mansinghka Vikash. Exploiting Conditional Independence for Efficient,
Automatic Multicore Inference for Church. — 2012.
[21] Blei David M, Ng Andrew Y, Jordan Michael I. Latent dirichlet allocation // the Journal
of machine Learning research. — 2003. — Vol. 3. — P. 993–1022.
[22] Wu Jeff. Reduced traces and JITing in Church : Ph. D. thesis / Jeff Wu ; Massachusetts
Institute of Technology. — 2013.
[23] CAPTCHA: Using hard AI problems for security / Luis Von Ahn, Manuel Blum,
Nicholas J Hopper, John Langford // Advances in Cryptology—EUROCRYPT 2003. —
Springer, 2003. — P. 294–311.
[24] Goodman N. D., Tenenbaum J. B. Probabilistic Models of Cognition. — 2014. — URL:
https://probmods.org/.

46

[25] A Repository for Generative Models (composite authors, edited by Andreas Stuhlmüller). —
2014. — URL: http://forestdb.org/.
[26] Wood Frank et al. Probabilistic Programming Tutorial (Machine Learning Summer School,
Iceland). — 2014. — URL: http://www.robots.ox.ac.uk/~fwood/anglican/teaching/
mlss2014/.
[27] Paige Brooks, Wood Frank. A Compilation Target for Probabilistic Programming
Languages. — 2014.
[28] Maddison Chris J, Tarlow Daniel. Structured Generative Models of Natural Source Code. —
2014.
[29] Liang Percy, Jordan Michael I, Klein Dan. Learning programs: A hierarchical Bayesian
approach. — 2010. — P. 639–646.
[30] Marsaglia George, Bray Thomas A. A convenient method for generating normal
variables. — 1964. — Vol. 6, no. 3. — P. 260–264.
[31] Box George EP, Muller Mervin E et al. A note on the generation of random normal
deviates. — 1958. — Vol. 29, no. 2. — P. 610–611.
[32] Devroye Luc. Non-Uniform Random Variate Generation. — 1986.
[33] Gulwani Sumit, Kitzelmann Emanuel, Schmid Ute. Approaches and Applications of
Inductive Programming (Dagstuhl Seminar 13502). — 2014. — Vol. 3, no. 12. — P. 43–
66. — URL: http://drops.dagstuhl.de/opus/volltexte/2014/4507.
[34] Looks Moshe. Program evolution for general intelligence // FRONTIERS IN ARTIFICIAL
INTELLIGENCE AND APPLICATIONS. — 2007. — Vol. 157. — P. 125.
[35] Muggleton Stephen, De Raedt Luc. Inductive logic programming: Theory and methods. —
1994. — Vol. 19. — P. 629–679.
[36] Muggleton Stephen, Feng Cao. Efficient induction of logic programs. — 1992. — Vol. 38. —
P. 281–298.
[37] De Raedt Luc, Kersting Kristian. Probabilistic inductive logic programming. — Springer,
2008.
[38] Kersting Kristian. An inductive logic programming approach to statistical relational

47

learning / IOS Press. — 2005. — P. 1–228. — URL: http://people.csail.mit.edu/
kersting/FAIAilpsrl/.
[39] Muggleton Stephen. Stochastic logic programs. — 1996. — Vol. 32. — P. 254–264.
[40] Exploiting compositionality to explore a large space of model structures / Roger Grosse,
Ruslan R Salakhutdinov, William T Freeman, Joshua B Tenenbaum. — 2012.
[41] Structure discovery in nonparametric regression through compositional kernel search /
David Duvenaud, James Robert Lloyd, Roger Grosse et al. — 2013.
[42] Hwang Irvin, Stuhlmüller Andreas, Goodman Noah D. Inducing probabilistic programs by
Bayesian program merging. — 2011.
[43] Gordon Andrew D. An Agenda for Probabilistic Programming: Usable, Portable,
and Ubiquitous // Workshop on Probabilistic Programming: Democratizing Machine
Learning. — 2013.
[44] Markov chain Monte Carlo without likelihoods / Paul Marjoram, John Molitor,
Vincent Plagnol, Simon Tavaré. — 2003. — Vol. 100, no. 26. — P. 15324–15328.
[45] Johnson Mark, Griffiths Thomas L, Goldwater Sharon. Adaptor grammars: A framework
for specifying compositional nonparametric Bayesian models. — 2007. — Vol. 19. — P. 641.
[46] Knuth Donald E. The Art of Computer Programming, 3rd edn., vol. 2. — 1998.
[47] Quinlan J. Ross. Simplifying decision trees. — 1987. — Vol. 27, no. 3. — P. 221–234.
[48] Bache K., Lichman M. UCI Machine Learning Repository. — 2013. — URL: http:
//archive.ics.uci.edu/ml.
[49] Tenenbaum Joshua B, Jonas Eric M, Mansinghka Vikash K. Stochastic Digital Circuits for
Probabilistic Inference. — 2008.
[50] Jonas Eric Michael. Stochastic architectures for probabilistic computation : Ph. D. thesis /
Eric Michael Jonas ; Massachusetts Institute of Technology. — 2014.

48

