Posledních několik let se hodně mluví a píše o datech a jejich analýze. Zejména o tom, jak rychle zpracovávat velká data pomocí moderních technologií postavených na MapReduce modelu. Existuje velké množství firem, které nabízejí služby zprostředkovávající vhled do dat pomocí vizualizací jednoduchých statistik. K dispozici je i řada knihoven a frameworků umožňujících vyvozovat z dat netriviální závěry s využitím standardních algoritmů strojového učení pro klasifikaci, shlukování apod. (za všechny jmenujme např. Apache Mahout). Relativně málo pozornosti je však věnováno obecným unsupervised nástrojům pro vyvozování informací z nestrukturovaných dat, kterých je většina. Přitom informace skryté v takových datech jsou často nejzajímavější.

Jedním z nejelegantnějších přístupů k unsupervised modelování dat jsou pravděpodobnostní programovací jazyky. Jedná se o vysokoúrovňové jazyky, které umožňují snadno popsat pravděpodobnostní model reprezentující data a automaticky z něho vyvozovat neznámé parametry. Interpret pravděpodobnostního programovacího jazyka se sám postará o inferenci, což programátorovi ušetří drahocenný čas strávený s odvozování vztahů a implementací kódu. Ve skutečnosti uživatel pravděpodobnostního jazyka vůbec nemusí umět programovat v běžném slova smyslu, ale vystačí si se znalostí statistiky a obecného charakteru svých dat.

Nejdůležitější částí kódu v pravděpodobnostních programovacích jazycích je návrh generativního modelu. Jedná se o pravděpodobnostní model, který v zásadě popisuje, jakým způsobem by mohla být vygenerována analyzovaná data, kdybychom znali všechny jeho parametry. Základními kameny takového modelu jsou vhodně poskládaná pravděpodobnostní rozdělení. Architektura modelu je volena tak, aby co nejlépe popisovala obecný charakter dat. Inferencí se rozumí proces nalezení takových parametrů, které maximalizují pravděpodobnost generování poskytnutých dat pomocí zvoleného modelu.

Na úvod začněme s jednoduchým ilustračním problémem. Mějme datový soubor, ve kterém jsou uložené výšky všech dospělých lidí v nějakém městě a řekněme, že nás zajímá průměrná hodnota μ a směrodatná odchylka σ těchto hodnot. Víme, že výška lidí je rozdělena podle normálního rozdělení

jak je ilustrováno na obrázku 1.

Proto stačí ve zvoleném jazyce zapsat, že je každá hodnota ve vstupních datech generována z normálního rozdělení a že chceme nalézt takové hodnoty jeho dvou parametrů, které jsou optimální vzhledem k analyzovaným datům. Tato úloha se dá samozřejmě řešit jednodušeji, neboť se vztahy pro výpočet těchto parametrů dají snadno odvodit. Existují však úlohy, kde je odvození obtížné, nebo dokonce není vůbec možné. Příkladem je následující rozšíření problému.

Mějme stejná vstupní data jako v předchozím příkladu, ale dejme si za cíl zjistit, jaká je průměrná výška mužů a průměrná výška žen v těchto datech. Na první pohled se může zdát, že se tyto údaje zjistit nedají, lze však využít obecných znalostí, které máme:

• jak výšky mužů, tak i výšky žen jsou rozděleny podle normálního rozdělení, ovšem s různými středními hodnotami
• muži jsou v průměru vyšší než ženy
• v populaci je v přibližně stejné množství mužů i žen

Jednou z možností, jak popsat data tak, abychom z nich byli schopní odvodit průměrnou výšku žen i mužů představuje následující generativní proces:

1. vyber hodnoty μ_muž, μ_žena tak, že μ_muž ≥ μ_žena
2. u každého člověka si hoď mincí a podle výsledku rozhodni, zda se jedná o muže nebo ženu
3. v případě, že jde o muže, vyber jeho výšku z normálního rozdělení N(μ_muž, σ²)
4. v případě, že jde o ženu, vyber její výšku z normálního rozdělení N(μ_žena, σ²)

Poté už stačí interpretu jazyka sdělit, aby nalezl takové hodnoty μ_muž, μ_žena a σ, které maximalizují pravděpodobnost dat při daném modelu.

Programovací jazyk JAGS

Jedním z nejpoužívanějších pravděpodobnostních programovacích jazyků v současnosti je jazyk JAGS. Jedná se o open source variantu jazyků WinBugs a OpenBugs, která používá stejnou syntax, ale je napsaná v jazyce C++. To umožňuje snadnou přenositelnost na různé platformy. Vzhledem k tomu, že JAGS nemá grafické uživatelské rozhraní a není v něm možné řešit předzpracování a následnou analýzu dat, je nejvhodnější ho použít jako knihovnu volanou z jiného programovacího jazyka. Nejčastější volbou je jazyk R, který je i syntakticky velmi podobný. Pro inferenci používá JAGS simulaci Markov Chain Monte Carlo (MCMC). Ta spočívá v tom, že jsou hodnoty neznámých náhodných proměnných po předchozí inicializaci postupně samplovány, čímž se jejich odhad neustále zpřesňuje. Po určitém počtu iterací algoritmu pravděpodobnostní rozdělení hodnot těchto proměnných zkonverguje. Poté je možné spustit další fázi algoritmu, ve které už se samplují jen proměnné, jejichž hodnoty chceme odvodit. Vzorky dat z této fáze se zprůměrují, čímž získáme očekávané hodnoty.

Pro použití v jazyce R je nejprve nutné nainstalovat interpret jazyka JAGS, který je k dispozici pro různé operační systémy, případně je možné jej zkompilovat přímo ze zdrojových kódů. Druhým krokem je instalace balíčku rjags do R. To je možné provést např. příkazem

install.packages("rjags")

přímo v interaktivním režimu interpretu jazyka R.

Nyní již máme připraveno vše k tomu, abychom mohli začít programovat. Vraťme se k prvnímu příkladu s výškami lidí ve městě. Vzhledem k tomu, že vzorek výšek lidí nemáme k dispozici, vytvořme si ho uměle:

N <- 10000
x <- rnorm(N, 173, 10)

Uvedený příklad vygeneruje náhodný vektor 10000 vzorků výšek z normálního rozdělení se střední hodnotou 173 a směrodatnou odchylkou 10. Naším cílem teď bude pouze na základě těchto náhodných dat zpětně odvodit neznámé parametry rozdělení. Kompletní kód v jazyce JAGS, který uložíme do souboru "vysky.bug", vypadá následovně:

model {
    for (i in 1:N) {
        x[i] ~ dnorm(mu, tau)
    }
    mu ~ dnorm(170, 0.0001)
    sigma ~ dunif(0, 50)
    tau <- pow(sigma, -2)
}

Každý kód v jazyce JAGS obvykle začíná klíčovým slovem model, které oznamuje, že bude následovat definice pravděpodobnostního modelu. Na řádcích 2-4 popisujeme generování výšek jednotlivých osob z normálního rozdělení s parametry mu a tau. Zdůrazněme, že v jazyce JAGS není druhým parametrem normálního rozdělení směrodatná odchylka, nýbrž převrácená hodnota rozptylu. Hodnoty mu jsou na řádku 5 generovány z normálního rozdělení. To je v tomto případě hrubým odhadem, který není informativní a má zanedbatelný vliv na výsledek inference. Podobně jsou hodnoty sigma generovány z rovnoměrného rozdělení z intervalu 0 až 50. To je opět neinformativní hrubý odhad, neboť skutečná směrodatná odchylka výšek lidí je určitě nižší než 50 cm. Konečně řádek 7 definuje vztah mezi tau a sigma. Všimněme si, že zatímco u generování z pravděpodobnostního rozdělení je použit symbol '~', u deterministického vztahu je to symbol '<-'.

Nyní už můžeme přistoupit k samotné MCMC simulaci v jazyce R:

library('rjags')

jags <- jags.model('vysky.bug', data = list('x' = x, 'N' = N))
update(jags, 1000)
jags.samples(jags, c('mu', 'sigma'), 1000)

Na prvním řádku nejprve načteme knihovnu 'rjags' do paměti. Na řádku 3 nahrajeme JAGS model ze souboru "vysky.bug" a přiřadíme vstupní data. Příkaz update na dalším řádku provede požadovaný počet iterací. Minimální počet iterací pro dosažení konvergence závisí jak na datech, tak na složitosti modelu a nedá se nijak jednoduše stanovit předem. U složitějších modelů je vhodné např. sledovat perplexitu na odložených datech a samplování zastavit ve chvíli, kdy perplexita přestane klesat. Pro tento jednoduchý příklad by však 1000 iterací mělo bez problému stačit. Konečně poslední příkaz samples provede dalších 1000 iterací, u kterých si však zaznamenává samplované hodnoty proměnných mu a sigma, které po skončení zprůměruje a vypíše na výstup. Tyto hodnoty by zhruba měly odpovídat parametrům, které byly použity pro umělé vygenerování vstupních dat, tedy mu ≈ 173 a sigma ≈ 10. Je důležité si uvědomit, že se jedná o stochastický proces, takže výsledek bude při každém spuštění trochu jiný.

Druhý příklad, který si klade za cíl odhadnout průměrné výšky mužů a žen je nepatrně složitější. Začněme vygenerováním umělých dat:

muzi <- rnorm(5000, 180, 10)
zeny <- rnorm(5000, 167, 10)
x <- sample(c(muzi, zeny))
N <- 10000

Nejprve vygenerujeme zvlášť vektor výšek mužů a zvláš'ť vektor výšek žen z normálních rozdělení. Obě rozdělení mají stejnou směrodatnou odchylku, ale různé střední hodnoty. Dále tyto seznamy spojíme a náhodně zamícháme, abychom si úlohu znesnadnili.

Nyní se podívejme na zápis modelu v jazyce JAGS:

model {
    for( i in 1 : N ) {
        T[i] ~ dcat(P[])
        x[i] ~ dnorm(mu[T[i]], tau)
    }
    theta ~ dnorm(0.0, 0.0001)
    mu[1] ~ dnorm(170.0, 0.0001)
    mu[2] <- mu[1] + abs(theta)
    P[1] <- 0.5
    P[2] <- 0.5
    sigma ~ dunif(0, 50)
    tau <- pow(sigma, -2)
}

Proměnná mu je tentokrát vektor délky 2, ve kterém první složka reprezentuje střední hodnotu výšky mužů a druhá střední hodnotu výšky žen. Pro každého člověka se nejprve vybere, zda se jedná o muže nebo o ženu. To je reprezentováno vektorem T délky N, kde 1 značí ženu a 2 muže. Tato rozhodnutí jsou náhodná, přičemž obě možnosti mají stejnou pravděpodobnost. Hodnota mu[1] je generována z neinformativního normálního rozdělení jako v předchozím příkladu, hodnota mu[2] je definována jako mu[1] + abs(theta), kde theta je z neinformativního rozdělení s nulovou střední hodnotou. Díky tomu, že je theta přičtena v absolutní hodnotě, je vždy dosaženo vztahu mu[1] ≤ mu[2]. U směrodatné odchylky výšek předpokládáme, že je stejná u mužů i u žen a je generována stejným způsobem jako dříve. Spuštění simulace je možné provést stejně jako v předchozím příkladu s tím rozdílem, že na výstupu v případě proměnné mu dostaneme místo skaláru vektor.

Na posledním příkladu si ukážeme, že se pravděpodobnostní programovací jazyky dají použít nejen pro přímočarou statistickou analýzu dat, ale i pro problémy, které se typicky řeší jinými způsoby. Takovou úlohou je například regrese. Cílem regresní analýzy je nalezení parametrů křivky tak, aby co nejlépe aproximovala nějakou množinu bodů. Křivkou může být například polynom třetího stupně

a hledanými parametry koeficienty w = (w₃, w₂, w₁, w₀), jak je ilustrováno na obrázku 2.

 

Modré body zde reprezentují prokládaná data, červená křivka je grafem hledaného polynomu. Typicky se tato úloha řeší hledáním nejmenší sumy čtverců chyb, lze na ni však nahlížet i Bayesovskou perspektivou. To, že modré body neleží přesně na červené křivce je obvykle způsobeno nějakým šumem, distribuovaným podle normálního rozdělení se střední hodnotou ležící na křivce. To je v obrázku znázorněno modrou křivkou. Lze tedy snadno navrhnou generativní model, který bude popisovat, jak by vznikla vstupní data, kdybychom znali parametry prokládané křivky.

Nejprve si opět vytvořme syntetická data:

x <- seq(-10, 10, by = 0.01)
N <- length(x)
w3 <- 1
w2 <- 0
w1 <- -1
w0 <- 10
epsilon <- rnorm(N, 0, 5)
y <- w3*x^3 + w2*x^2 + w1*x + w0 + epsilon

Z kódu je vidět, že jsou body rozprostřeny podle normálního rozdělení se směrodatnou odchylkou 5 kolem křivky s rovnicí

Následuje popis generativního modelu v jazyce JAGS, který uložme do souboru 'regression.bug':

model {
    for (i in 1:N) {
        y[i] ~ dnorm(y0[i], tau)
        y0[i] <- w3*pow(x[i], 3) + w2*pow(x[i], 2) + w1*x[i] + w0
    }
    w3 ~ dnorm(0, 0.0001)
    w2 ~ dnorm(0, 0.0001)
    w1 ~ dnorm(0, 0.0001)
    w0 ~ dnorm(0, 0.0001)
    tau <- pow(sigma, -2)
    sigma ~ dunif(0, 100)
}

Vidíme, že jsou datové body generovány z normálního rozdělení se střední hodnotou w3*pow(x[i], 3) + w2*pow(x[i], 2) + w1*x[i] + w0 a směrodatnou odchylkou sigma. Parametry w3, w2, w1, w0 jsou generovány z neinformativního normálního rozdělení. Pro úplnost doplňme kód pro nalezení optimálních hodnot těchto parametrů:

jags <- jags.model('regression.bug', data = list('x' = x, 'y' = y, 'N' = N))
update(jags, 1000)
jags.samples(jags, c('w3', 'w2', 'w1', 'w0'), 1000)

Cílem představení uvedených příkladů bylo ilustrovat způsob použití pravděpodobnostních programovacích jazyků. Ve skutečnosti existují pro všechny tři úlohy efektivnější řešení. Ta však vyžadují odvozování vztahů a výrazně větší množství programového kódu. Některé modely jako např. LDA navíc efektivněji řešit nejdou (nebo alespoň žádné takové řešení není známé) a pravděpodobnostní programovací jazyky jsou tak nejlepší volbou pro jednoduché experimenty. Pro další inspiraci lze použít jako zdroj příkladů archiv na webu projektu OpenBugs.

Zdroje obrázků:

• http://research.microsoft.com/en-us/um/people/cmbishop/PRML/webfigs.htm
• http://www.sentientdevelopments.com/2006/01/economist-bayes-helps-explain-how-mind.html