Markov-lánc Monte Carlo

Ahogy a Bayesian szemlélettel foglalkozó bejegyzésben említettem, végtelen nagyságú tér esetén mintavételt kell alkalmaznunk ahhoz, hogy a teret leírjuk. A mai bejegyzés egy olyan mintavételi eljárásról szól, a Markov-lánc Monte Carloról (Markov chain Monte Carlo, MCMC ), amit gyakran alkalmaznak erre a feladatra, és nem csak Bayesian hanem más területen is pl: megerősítéses tanulás.


Kezdjük talán azzal, miért terjedt el ez a mintavételi mód? A legegyszerűbb mintavételi eljárás az úgynevezett Hálós mintavétel1. Ez szabályos, általában egyforma nagyságú cellákra osztja fel a teret és mindegyik cellából ugyanannyi mintát vesz. Ha csak a tér egy részét szeretnénk megismerni, ez a módszer elég erőforrás pazarló. A Bayesian analízis során általában ez a helyzet. Az átlagtól távolodva egyre kisebb egy adott Hipotézis valószínűsége, és egy idő után ez a valószínűség értelmezhetetlenül kicsi lesz. Nem optimális ezekről a helyekről is mintát venni. Inkább azt szeretnénk, hogy az értelmezhető tartományon belüli területet jobban részletezzük. A Markov-lánc Monte Carlo pont ezt csinálja. A Hipotézistér pontjait annál gyakrabban vizsgálja meg, minél nagyobb a valószínűsége, hogy az adott Hipotézis igaz.

Az MCMC két módszer kombinációja: a Markov-lánc és a Monte-Carlo-módszer. Nézzük mik ezek.

Monte-Carlo-módszer

Kezdjük egy kis érdekességek: miért Monte-Carlo a neve? A Monte-Carlo-módszert az 1940-es években fejlesztette ki Stanislaw Ulam. Rögtön a publikálás után Neumann János felismerte a módszer fontosságát, és az első programozható számítógépre, az ENIAC-ra meg is írta a maga változatát. Viszont Neumann és Ulam közös munkája titkos volt, így kódnevet adtak neki, ez lett a Monte-Carlo. Azért választották ezt, mert Ulam sok pénzt vesztett itt szerencsejátékon, és mint később látni fogjuk a véletlenszerűség szerves része ennek a módszernek. Miért volt titkos ez a kutatás? A Monte-Carlo-módszer a Manhattan terv szimulációs részéhez kapcsolódót.2

Ulam alapötlete az volt, hogy ha egy számítást nem tudunk megoldani analitikusan, akkor csak szimuláljuk és végül számoljuk össze az eredményt. Pl: nem tudjuk a Pi-t kiszámolni. Akkor rajzoljunk egy kört egy négyzetbe majd dobjunk véletlenszerűen pontokat ebbe a négyzetbe. Ha megszámoljuk mennyi van a körön belül már közelítőleg tudjuk a Pi értékét. Ha növelni akarjuk a pontosságot, csak használjunk több pontot.

A módszer mai szemmel triviálisnak tűnik, de ne felejtsük el, hogy Ulam elmélete akkor született amikor még nem voltak programozható számítógépek. Ennek megfelelően az általunk használt véletlenszám generáló eljárás se létezett. Pont Neumann érdeme, hogy először sikerült egy olyan algoritmust kidogozni, ami már értékelhető időn belül volt képes véletlenszerű számot generálni.

A módszer általános lépései:

  1. Meghatározunk egy domént a lehetséges bemenetre
  2. Véletlenszerű mintát veszünk ebből a doménből valamiféle sűrűségfüggvény felhasználásával
  3. Elvégzünk egy determinisztikus számítást a 2. pont eredményén
  4. Összegezzük az eredményt

Nézzünk egy példát. Szeretnénk tudni mekkora a valószínűsége, hogy -1,96 és +1,96 között van egy standard normál eloszlás értéke:

import numpy as np
# mintavetelek szama
runs = 100000
# 1. lépés -- domain
atlag = 0
std = 1
# 2. lépés -- mintavetel
simulation = np.random.normal(atlag,std,runs)
# 3. lépés -- determinisztikus számítás
szamitas= [ x for x in simulation if x >= -1.96 and x <= 1.96 ]
# 4. lépés -- összegzés
osszegzes=  len(szamitas)/runs
print(osszegzes)

Mint tudjuk a fenti eredménye kb. 95% a valóságban, és ahogy vártuk a fenti kód is hasonló számot3 ad.

Most nézzük a Markov-láncot.

Markov-lánc

A Markov láncok elnevezésének nincs olyan érdekes története mint a Monte-Carlo-módszernek. A nevét a kidolgozója után, Andrej Andrejevics Markov orosz matematikus után kapta. Maga az eljárás külön bejegyzést érdemel, amit nemsokára tervezek, így itt nem fogom részletezni az általános elvét. Most egyelőre elég annyit tudni róla, hogy egy Markov-lánc esetén a jövőbeli állapot csak a jelenlegi állapottól függ. Az nem számít, hogy jutottunk el ehhez a jelenlegi állapothoz. Vagyis a Markov-lánc memória nélküli.

Metropolis-Hasting algoritmus

A Markov-lánc általános elmélete helyett nézzünk inkább a MCMC egy gyakorlati megvalósítását, a Metropolis-Hasting algoritmust. Miért pont ezt? Mert ez az az algoritmus, amit az általunk használt Bayesian szoftverek, a STAN és a PyMC3, használnak.

Az algoritmust Nicholas Metropolis után lett elnevezve, aki 1953-ban publikálta a Equation of State Calculations by Fast Computing Machines cikkét. A cikk társzerzői többek között Teller Ede és szintén magyar felesége Teller Auguszta Mária. Mária írta mellesleg a MCMC első kódját is.

Az algoritmus a következő lépéseket hajtja végre:

  1. Választunk egy kezdő értéket x_i számára. Ez lehet véletlenszerű, de nem muszáj.
  2. Véletlenszerű mintát veszünk egy egyszerű eloszlásból, például uniform eloszlásból, vagy Gaussianból. Ez a Monte-Carlo része az algoritmusnak. Ez a minta lesz az x_{i+1} mintánk, ha elfogadjuk a következő két lépésben.
  3. Kiszámítjuk mi a a valószínűsége annak, hogy elfogadjuk a 2. pontban vett mintát. Ez a Metropolis-Hasting kritériummal történik. Hogy mi is ez lsd. lentebb.
  4. Ha a 3. pontban számított szám nagyobb mint egy véletlenszerű minta az uniform[0,1] eloszlásból, akkor elfogadjuk a 2. pontban vett mintát, ha nem akkor dobjuk és helyettesítjük a x_i   -vel. Vagyis ebben az esetben a x_{i+1} = x_i    . A 3. és ez a pont a Markov-lánc része az algoritmusnak.
  5. Addig ismételjük a 2. pont és 4. pont közötti részt amíg elegendő mintát gyűjtünk

Ez nagyjából elég egyértelmű, a 3. pontban említett Metropolis-Hasting kritérium kivételével.4 Nézzük mi is ez:

(1)   A = p_a(x_{i+1}|x_i) = min\left(  1, \frac{p(x_{i+1}) \cdot q(x_i|x_{i+1}) }{ p(x_{i}) \cdot q(x_{i+1}|x_{i}) } \right)

Ahol:

  • p(x_{i}) — a x_i   minta előfordulásának valószínűsége a kereset modellben
  • q(x_i|x_{i+1}) — a x_i   minta előfordulásának valószínűsége a ha a x_{i+1}   pontban állunk. Magyarul menyire elképzelhető, hogy a x_{i+1}    pontból a x_i   pontba lépünk. Értelemszerűen ez attól függ, hogy milyen eloszlást használunk a következő pont kiválasztására. Ha Gaussiant akkor a közeli pontoknak nagyobb az esélye. Ha Uniformot akkor mindegyiknek egyforma az esélye, stb.
  • p(x_{i+1}) — a x_{i+1}   minta előfordulásának valószínűsége a kereset modellben
  • q(x_{i+1}| x_i ) -- a latex x_{i+1} &s=2 $ minta előfordulásának valószínűsége a ha a x_{i}   pontban állunk.

Ha megnézzük a fenti kifejezést, akkor lényegében azt mondja, hogy ha az új pont valószínűsége nagyobb, mint a korábbi ponté, akkor mindig elfogadjuk. Ha kisebb akkor a kiválasztás valószínűsége aranyos a pontok valószínűségével. Miért jó ez? Miért nem dobjuk egyszerűen a kisebb pontokat? A válasz: ezzel az eljárással el tudjuk kerülni, hogy lokális minimumokba ragadjunk.

Ha a q(x_{i+1}| x_i ) =   q(x_i|x_{i+1})  , vagyis annak az esélye, hogy x_i   közötti mozgás x_{i+1}   egyformán valószínű, akkor a (1) egyszerűsíthetjük:

(2)   A = min\left(  1,  \frac{p(x_{i+1}) }{ p(x_{i})   } \right)

Nézzünk egy naiv Python megvalósítást ahol az adatok egy \lambda = 1 exponenciális eloszlásból származnak. Normál eloszlást fogunk használni a 2. lépésben, így használhatjuk a (2)-t:

# MCMC
import numpy as np
def keresett_eloszlas(x):
    if x  <= 0:
        return 0
    else:
        return np.exp(-x)
n = 10000
# 1. lépés -- iniciálás
x = [3]
for i in range(n):
    jelenlegi_x = x[-1]
    # 2. lépés -- minta egy standard normál eloszlásból
    felteles_x = jelenlegi_x + np.random.normal(loc=0.0, scale=1.0, size=1)[0]
    # 3.lépés --  Metropolis-Hasting kritérium számítása
    try:
        A = keresett_eloszlas(felteles_x)/keresett_eloszlas(jelenlegi_x)
    except:
        continue
    # 4. lépés -- minta egy uniform eloszlásból és az A összehasonlítása ezzel
    uniform_minta = np.random.uniform(0,1,1)[0]
    if uniform_minta < A:
        x.append(felteles_x)
    else:
        x.append(jelenlegi_x)
# Ábrázoljuk
import matplotlib.pyplot as plt
bins_tavolsag = 10
bins = [  x/bins_tavolsag  for x in range(8*bins_tavolsag) ][1:]
weights = np.ones_like(x)/float(len(x))
plt.hist(x, bins = bins , density=True, label="MCMC minta")
plt.plot( bins, [ keresett_eloszlas(bin) for bin in bins ], label="Valós eloszlás" )
plt.title("Markov-lánc Monte Carlo példa")
plt.xlabel("Megfigyelés")
plt.ylabel("Valószínűség")
plt.legend()
plt.show()

Aminek az eredménye:

Mint látható a MCMC mintavétel elég jól közelíti a valós eloszlást.

Lábjegyzet

  1. angolul: grid sampling
  2. Bővebben: Computing and the Manhattan Project
  3. Nekem konkrétan 94,915% lett az eredmény, de mivel véletlenszerű folyamatról beszélünk ez minden egyes próbával egy kicsit változhat.
  4. Angolul “acceptance probabilty”-nek is nevezik néha

Irodalom

Hírdetés

Markov-lánc Monte Carlo” bejegyzéshez 2 hozzászólás

Vélemény, hozzászólás?

Adatok megadása vagy bejelentkezés valamelyik ikonnal:

WordPress.com Logo

Hozzászólhat a WordPress.com felhasználói fiók használatával. Kilépés /  Módosítás )

Twitter kép

Hozzászólhat a Twitter felhasználói fiók használatával. Kilépés /  Módosítás )

Facebook kép

Hozzászólhat a Facebook felhasználói fiók használatával. Kilépés /  Módosítás )

Kapcsolódás: %s