Імовірнісні моделі: семплірованіє

    І знову здрастуйте! Сьогодні я продовжую серію статей в блозі Surfingbird, присвячену різним методам рекомендацій, а також іноді і просто різного роду імовірнісним моделям. Давним-давно, здається, минулої п'ятниці влітку минулого року, я написав невеликий цикл про графічних імовірнісних моделях: перша частина вводила основи графічних імовірнісних моделей, під другій частині було кілька прикладів, частина 3 розповідала про алгоритм передачі повідомлень, а в четвертої частини ми коротко поговорили про варіаційних наближеннях. Цикл закінчувався обіцянкою поговорити про Семплірування — ну що ж, не минуло й року. Взагалі кажучи, в цьому міні-циклі я поведу мову більш предметно про модель LDA і про те, як вона допомагає нам робити рекомендації текстового контенту. Але сьогодні почну з того, що виконаю давню обіцянку і розповім про Семплірування в імовірнісних моделях — одному з основних методів наближеного висновку.
 
 
 
 У чому завдання
Насамперед я нагадаю, про що ми взагалі говоримо. Один з основних інструментів машинного навчання — це теорема Байеса:
 
Тут D — це дані, θ — параметри моделі, які ми хочемо навчити, а — це розподіл θ за умови наявних даних D ; про теорему Байєса ми вже докладно говорили одній з попередніх статей . У машинному навчанні ми, як правило, хочемо знайти апостеріорне розподіл , а потім передбачити нові значення змінних . У складних графічних моделях все це, як ми обговорювали в попередніх серіях, зазвичай зводиться до того, що у нас є велике і заплутане розподіл різноманітних випадкових величин , яке розкладається в добуток розподілів попроще, і наше завдання — провести маргіналізацію , т. е. підсумовувати по частині змінних або знайти сподівання функції від частини змінних. Зауважимо, що всі наші завдання так чи інакше зводяться до того, щоб підрахувати математичне сподівання різних функцій за умови складного розподілу, зазвичай умовного розподілу моделі, в яку підставлені значення деяких змінних. Також можна вважати, що ми вміємо рахувати значення спільного розподілу в будь-якій його точці — це зазвичай нескладно випливає з загального вигляду моделі, а складність саме в тому, щоб по купі цих змінних підсумувати-проинтегрировать.
 
Ми вже знаємо, що точний висновок вести складно, та ефективні алгоритми виходять тільки для випадку фактор-графа без циклів або з маленькими циклами. Проте реальні моделі часто містять купу змінних, які пов'язані один з одним досить щільно і утворюють масу циклів. Ми трохи говорили про варіаційних наближеннях — одному можливий спосіб подолати виникаючі труднощі. А сьогодні я розповім про інше, ще більш популярному методі — про Семплірування.
 
 Наближення Монте-Карло
Ідея проста до надзвичайності. Ми вже з'ясували, що в принципі все, що нам потрібно, виражається у вигляді очікувань різноманітних функцій по нашому складному і заплутаному розподілу p (x ). Як підрахувати сподівання функції з розподілу? Якщо ми вміємо семпліровалісь з цього розподілу, тобто брати випадкові точки з розподілу p (x ), то очікування будь-якої функції можна наблизити як середнє арифметичне в семплірованних точках:
 
Наприклад, щоб порахувати середню (матожіданіе), потрібно усереднити по семплам x .
 
Таким чином, вся наша задача фактично зводиться до того, щоб навчитися семпліровалісь точки з розподілу за умови, що ми можемо в будь-якій точці вважати значення цього розподілу.
 
 У чому тут складність
Людині, вихованому на функціях від однієї змінної (тобто фактично будь-якій людині, прослухають базовий курс матаналізу у вузі), може здатися, що великої проблеми тут немає, і всі ці розмови ні про що — відомо ж, як наближено порахувати інтеграл, треба просто розбити відрізок на частини, порахувати функцію в середині кожного маленького отрезочка, прямокутники там, трапеції, ну ви розумієте… І дійсно, в Семплірування з розподілів від однієї змінної немає нічого складного. Складнощі, як завжди, починаються з ростом розмірності — якщо в розмірності 1 потрібно порахувати 10 значень функції, щоб розбити інтервал на відрізки зі стороною 0.1, то в розмірності 100 таких значень потрібно буде вже 10 100 , що набагато сумніше. Цей та інші ефекти, які роблять роботу в високої розмірності абсолютно несхожою на звичну нам розмірність 2-3, отримали назву «прокляття розмірності» (curse of dimensionality); там є й інші цікаві контрінтуітівное ефекти, можливо, варто коли-небудь поговорити про це докладніше.
 
Є і ще одна складність — вище я писав, що ми можемо вважати значення розподілу в будь-якій точці. Насправді ж зазвичай ми можемо рахувати не значення розподілу p (x ), а значення деякої функції p * (x ), яка пропорційна p (x ). Згадайте теорему Байєса — з неї ми отримуємо не саме апостеріорне розподіл, а пропорційну йому функцію:
 
Зазвичай на цьому місці ми говорили, що цього достатньо, а точне значення ймовірності можна просто підрахувати, нормалізувавши отриманий розподіл так, щоб воно підсумовані в одиницю (тобто обчисливши знаменник в теоремі Байеса). Однак це якраз і є та складне завдання, яку ми зараз намагаємося навчитися вирішувати: підсумувати всі значення складного розподілу. Тому в реальному житті ми не можемо розраховувати на значення істинних ймовірностей, тільки на деяку пропорційну їм функцію p * .
 
 Що ж робити: семпліруем під графіком
Незважаючи на всі перераховані вище складності, семплірованіє робити все-таки можливо, і це дійсно один з головних прийомів наближеного висновку. Більш того, як ми скоро побачимо, робити семплірованіє зовсім не так складно, як здається. У цьому розділі я розповім про основну ідею, на якій так чи інакше засновані всі алгоритми семплювання. Ідея така: якщо рівномірно вибрати точку під графіком функції щільності розподілу p , то її X-координата буде взята саме з розподілу p . Це дуже легко зрозуміти інтуїтивно: уявіть собі графік щільності і розбийте його на маленькі стовпчики, приблизно так (графіки зроблені в matplotlib):
 
Тепер видно, що коли ви берете випадкову точку під графіком, кожна X-координата буде зустрічатися з імовірністю, пропорційною висоті свого стовпчика, тобто якраз із значенням функції щільності на цьому стовпчику.
 
Таким чином, нам потрібно тільки навчитися брати випадкову точку під графіком щільності потрібного розподілу. Як же це зробити?
 
Перша ідея, яка реально працює і часто застосовується на практиці — так звана вибірка з відхиленням (rejection sampling ). Давайте припустимо, що у нас є якесь інший розподіл q (точніше, пропорційна йому функція q * ), яке ми вже вміємо семпліровалісь. Зазвичай це небудь класичний розподіл — рівномірний або нормальне. Тут варто відзначити, що семпліровалісь з нормального розподілу — це теж не цілком тривіальна задача, але ми зараз її пропустимо, для нас зараз досить знати, що «люди так робити вміють». Далі, припустимо, що ми знаємо про нього таку константу c , що для всіх x . Тоді ми зуміємо семпліровалісь p таким чином:
 
     
беремо семпл x з розподілу q * (x );
 беремо випадкове число u рівномірно з інтервалу ;
 обчислюємо p * (x ); якщо воно більше u , то x додається в семпли, а якщо менше (тобто якщо u не потрапило під графік щільності p * ), то x відхиляється (звідси і вибірка з відхиленням).
 
 
Ось картинка, що ілюструє вибірку з відхиленням (нормальний розподіл, підібране так, щоб мажоріровать суміш двох інших гауссіанов):
 
Якщо точка під графіком cq * потрапляє під графік p * , тобто в синю зону, ми її беремо; а якщо над ним, тобто в червону зону, — не беремо. Знову ж, легко інтуїтивно зрозуміти, чому це працює: ми беремо випадкову точку рівномірно під графіком cq * , а потім залишаємо тільки ті точки, які при цьому потрапили під графік p * — зрозуміло, що при цьому у нас залишиться якраз рівномірний розподіл під графіком p * . Настільки ж зрозуміло, що якщо нам вдасться вибрати q досить схожим на p (наприклад, ми знаємо приблизно, де зосереджено p , і накриваємо цю область гауссіаном), то цей метод буде набагато ефективніше, ніж просто розбивати весь простір на однакові шматочки. Але для того, щоб метод працював, потрібно, щоб cq дійсно досить добре наближало p — якщо площа синьої зони становить 10 -10 від площі червоної, ніяких семплів не вийде…
 
Є різні варіанти цієї ідеї, на яких я не буду зупинятися детально. Якщо коротко, то ця ідея і її варіації добре працюють зі складними плотностями в розмірності, що обчислюється одиницями, але в дійсно високої розмірності таки починаються складності. Через ефектів прокляття розмірності класичні пробні розподілу q отримують різні неприємні властивості (наприклад, в високої розмірності шкірка апельсина займає майже весь його обсяг, і, зокрема, нормальний розподіл теж майже цілком зосереджено в дуже тонкій «шкірці», а зовсім не біля нуля), семпли в потрібних місцях не вдається отримати, і вся справа часто провалюється. Потрібні інші методи.
 
 Алгоритм Метрополіса-Гастингса
Суть цього алгоритму заснована на тій же ідеї: ми хочемо взяти точку рівномірно під графіком функції. Однак підхід тепер інший: замість того щоб намагатися накрити весь розподіл «ковпаком» функції q , ми будемо діяти зсередини: будемо будувати випадкове блукання під графіком функції, переходячи від однієї точки до іншої, і час від часу брати поточну точку блукання в якості семплу. Оскільки в даному випадку це випадкове блукання є марковской ланцюгом (тобто його наступна точка залежить тільки від попередньої, а пам'яті ніякої нету), цей клас алгоритмів називається ще MCMC-семплюванням (Markov chain Monte Carlo). Щоб довести, що алгоритм Метрополіса-Гастингса дійсно працює, потрібно знати властивості марковських ланцюгів; але ми зараз нічого формально доводити не збираємося, так що я просто цнотливо залишу посилання на вікі .
 
Фактично ж алгоритм схожий на все ту ж вибірку з відхиленням, але з важливою відмінністю: раніше розподіл q було єдине для всього простору, а тепер воно буде локальним і буде змінюватися з часом, залежатиме від поточної точки блукання. Як і колись, ми розглядаємо сімейство пробних розподілів , де — поточний стан. Але тепер q не повинно бути наближенням p , а повинно просто бути яким-небудь семпліруемим розподілом, яке зосереджено в околиці поточної точки — наприклад, добре підходить сферичний Гауссіан з невеликою дисперсією. Кандидат в новий стан x 'тепер буде семпліровалісь з . Чергова ітерація алгоритму починається з стану і полягає в наступному:
 
     
вибрати x 'з розподілу ;
 обчислити
 
(Не треба лякатися, для симетричних розподілів, наприклад гауссіана, дорівнює 1, це просто поправка на асиметрію);
 З імовірністю a (1, якщо a > 1) , інакше .
 
 
Суть в тому, що ми переходимо в новий центр розподілу, якщо приймемо черговий крок, і виходить випадкове блукання, залежне від розподілу p : ми завжди беремо крок, якщо в цей бік функція щільності збільшується, і іноді відкидаємо, якщо зменшується (відкидаємо не завжди, тому що нам треба вміти виходити з локальних максимумів і блукати під усім графіком p ). Зауважу, до речі, що коли крок відкидається, ми не просто відкидаємо нову точку, а повторюємо другий раз те ж саме x (t ) — таким чином висока ймовірність повторювати семпли навколо локальних максимумів розподілу p , що відповідає більшої щільності.
 
І знову картинка — ось типове блукання на площині під графіком суміші двох двовимірних гауссіанов:
 
Щоб підкреслити, наскільки це насправді простий алгоритм, я приведу код, яким був породжений цей графік.
 Код на python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.path as path
import matplotlib.mlab as mlab
import scipy.stats as stats

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

z1 = stats.multivariate_normal([0,0],[[1.0,0],[0,1.0]])
z2 = stats.multivariate_normal([-2,-2],[[1.5,0],[0,0.5]])
def z(x): return 0.4*z1.pdf(x) + 0.6*z2.pdf(x)

Z1 = mlab.bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = mlab.bivariate_normal(X, Y, 1.5, 0.5, -2, -2)
Z = 0.4*Z1 + 0.6*Z2

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

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

p = path.Path(samples,codes)

fig, ax = plt.subplots()
ax.contour(X, Y, Z)
ax.add_patch(patches.PathPatch(p, facecolor='none', lw=0.5, edgecolor='gray'))
plt.show()

 
Отже, ми отримали випадкове блукання під графіком щільності p (x ); насправді це ще треба довести, але в математику ми заглиблюватися не будемо, так що повірте на слово. Що тепер робити з цим випадковим блуканням, нам же потрібні незалежні семпли, а тут виходить, що наступна точка блукання завідомо лежить поруч з попередньою, і ніякої незалежності немає і в помині?
 
Відповідь проста: потрібно брати не кожну точку, а тільки деякі, відокремлені один від одного багатьма кроками. Оскільки це випадкове блукання, то якщо більша частина q зосереджена в радіусі ε, а загальний радіус, в якому зосереджено p , дорівнює D , то для отримання незалежного семплу потрібно буде порядку кроків. У це ви, знову ж, можете повірити на слово, а можете згадати класичну задачку з курсу теорії ймовірностей про те, як далеко за n кроків піде точка, яка з однаковою ймовірністю рухається наліво і направо; піде вона в середньому на , так що логічно, що кроків буде квадратичне число.
 
А тепер головна фішка: це квадратичне число залежить від радіусів двох розподілів, але при цьому зовсім не залежить від розмірності . Семплірування за методом марковських ланцюгів Монте-Карло в набагато меншому ступені схильне прокляттю розмірності, ніж інші методи, які ми обговорювали; саме тому воно і стало, можна сказати, золотим стандартом. На практиці, правда, D оцінити важко, так що в конкретних реалізаціях зазвичай роблять так: спочатку деякий досить число перших кроків, скажімо 1000, відкидають (це називається burn-in, і це дійсно важливо, тому що початкова точка може потрапити в невдалу область p , так що перед початком процесу треба точно переконатися, що ми вже ходили досить довго), а потім беруть кожен n -й семпл, де n підбирають експериментально, виходячи з реально получающейся автокореляції в послідовності семплів.
 
 Семплірування по Гиббсу
І останній на сьогодні алгоритм семплювання — семплірованіє по Гиббсу. Це насправді окремий випадок алгоритму Метрополіса-Гастингса, дуже популярний і простий окремий випадок.
 
Ідея семплювання по Гиббсу ну зовсім проста: припустимо, що ми знаходимося в дуже великої розмірності, вектор x дуже довгий, і нам складно вибирати весь семпл відразу, не виходить. Давайте спробуємо вибирати семпл не весь відразу, а покомпонентно. Тоді напевно ці одномірні розподілу виявляться простіше, і семпл ми оберемо.
 
Формально кажучи, ми пробігаємо по компонентах вектора , на кожному кроці фіксуємо всі змінні, крім однієї, і вибираємо послідовно з розподілу
 
 
Від загального випадку блукання по Метрополіс-Гастінгс це відрізняється тим, що тепер ми рухаємося на кожному кроці вздовж однієї з координатних осей; ось відповідна картинка:
 
 
Семплірування по Гиббсу — це окремий випадок алгоритму Метрополіса для розподілів , і ймовірність прийняття кожного семплу виходить завжди дорівнює 1 (можете це довести в якості вправи). Тому семплірованіє по Гиббсу сходиться, і, так як це таке ж випадкове блукання по суті, вірна та ж квадратична оцінка.
 
Для семплювання по Гиббсу не потрібно ніяких особливих припущень або знань. Можна швидко зробити працюючу модель, і тому це дуже популярний алгоритм. Зауважимо ще, що у великих розмірностях може виявитися ефективніше семпли по кілька змінних відразу, а не по одній — наприклад, часто буває, що у нас двочастковий граф з змінних, в яких всі змінні з однієї частки пов'язані з усіма змінними з іншої частки (ну або з багатьма), а між собою не пов'язані. У такій ситуації сам Бог велів зафіксувати всі змінні однієї частки і просемпліровать всі змінні в іншій частці одночасно (це можна розуміти буквально — оскільки при такій структурі всі змінні однієї частки умовно незалежні за умови іншого, їх можна семпліровалісь незалежно і паралельно), потім зафіксувати всі змінні другого частки і так далі.

Висновок
Сьогодні ми встигли досить багато: ми поговорили про різні методи семплювання, яке є одним з основних інструментів наближеного виводу в складних імовірнісних моделях. Далі я планую розвивати цей новий міні-цикл в сторону тематичного моделювання, точніше, моделі LDA (latent Dirichlet allocation, латентне розміщення Дирихле). У наступній серії я розповім про саму модель (коротко я її вже описував , тепер можна поговорити більш детально) і поясню, як для LDA працює семплірованіє по Гиббсу. А потім ми будемо поступово переходити до того, як можна застосовувати LDA і її численні розширення для рекомендаційних систем.

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.