Bayesovský bandita: chytřejší a levnější A/B testování

Slot Machines

Představme si jednoduchou situaci. Řekněme, že jsme vlastníky internetového e-shopu s radiátory a chceme začít prodávat nejnovější typ. Ten má sice o něco lepší vlastnosti než staré modely, ale je o třetinu dražší, a zákazníci tak stále preferují starší varianty. Proto se rozhodneme, že ho začneme propagovat na úvodní stránce prostřednictvím banneru. Na prodeji tohoto modelu nám opravdu hodně záleží, a tak jsme si zaplatili profesionála, který připravil tři grafické návrhy. Otázka teď zní: jaký grafický návrh vybereme?

Nejjednodušší možností, kterou v takovém případě volí většina majitelů e-shopů, je vlastní úsudek. Ten je ovšem velmi subjektivní a vůbec nemusí odpovídat preferencím zákazníků. Lepší možností je využití takzvaného A/B testování. To spočívá v tom, že po nějakou dobu všechny bannery střídáme a pro každý z nich měříme poměr prokliků vůči počtu jeho zobrazení (CTR). Po nasbírání průkazného množství dat vybere banner s nejvyšší hodnotou CTR.

Tato metoda však není ideální. Návštěvnost e-shopu s radiátory asi nebude nijak závratná, proto může trvat velmi dlouho než nasbíráme dostatek dat pro kvalifikované rozhodnutí. Pokud bude navíc některý z bannerů výrazně méně atraktivní pro uživatele, jeho zobrazováním namísto jiného můžeme přijít o velké peníze. Naším cílem je tedy odhalit nejatraktivnější banner co nejdříve.

Moje oblíbená metoda, která patří mezi nejefektivnější, nese název bayesovský bandita. Pojmenování pochází z optimalizační úlohy známé jako problém mnohorukého bandity. Jedná se o problém, ve kterém máme řadu výherních automatů (v angličtině nazývaných one-armed bandit), každý s jinou pravděpodobností výhry, a naším cílem je nalézt posloupnost, v jaké na jednotlivých automatech hrát tak, abychom maximalizovali výnos. Jinými slovy, chceme co nejrychleji nalézt automat, u kterého je pravděpodobnost výhry nejvyšší a dále hrát pouze na tomto automatu. Analogie s bannery v e-shopu je přímočará − naším cílem je co nejdříve nalézt banner s nejvyšším CTR.

Jak sám název metody napovídá, je k tomu využita bayesovská statistika. Konkrétně jsou CTR jednotlivých bannerů odhadována samplováním z beta rozdělení. V obecnosti se jedná o rozdělení s parametry α a β definované jako

Beta(x|\alpha, \beta) = \dfrac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1 - x)^{\beta-1}.

Pro ilustraci jsou křivky hustoty pravděpodobnosti beta rozdělení pro některé hodnoty α a β jsou zachyceny na obrázku 1.

Pro porozumění je dobré si uvědomit, že beta rozdělení je v podstatě pravděpodobnostní rozdělení popisující pravděpodobnosti pravděpodobností. V našem případě je to rozdělení přes pravděpodobnosti prokliků (CTR). Tvar křivky určující beta rozdělení je dán parametry α a β. Pokud jsou oba rovny 1 (černá křivka), všechny hodnoty CTR jsou stejně pravděpodobné. V případě volby odlišných hodnot parametrů dochází k tomu, že některé hodnoty CTR jsou více pravděpodobné než jiné. Obecně platí, že je střední hodnota beta rozdělení určena vztahem

E[x] = \dfrac{\alpha}{\alpha + \beta}.

Pokud tedy budou například α i β stejné a zároveň větší než 1, bude vždy nejpravděpodobnější hodnotou CTR číslo kolem 0.5. Pokud budou hodnoty α a β v poměru 25:75, bude střední hodnota CTR 1/4 (modrá křivka). Další zajímavou vlastností beta rozdělení je fakt, že čím jsou hodnoty parametrů vyšší, tím menší je jeho rozptyl. Z praktického hlediska to znamená, že čím vyšší hodnoty parametrů α a β ve správném poměru zvolíme, tím přesnější odhad CTR budeme mít (zelená a červená křivka).

Nyní se můžeme pustit do popisu algoritmu bayesovského bandity. Hlavní myšlenka spočívá v tom, že CTR každého banneru modelujeme pomocí beta rozdělení. Při hledání vhodného banneru k zobrazení náhodně vylosujeme hodnotu CTR ze všech tří beta rozdělení a zobrazíme ten banner, jehož vylosovaná hodnota CTR bude nejvyšší. Pokud uživatel na banner klikne, zvýšíme parametr α odpovídajícího beta rozdělení o 1. Pokud uživatel neklikne, zvýšíme o 1 hodnotu parametru β. Tento postup opakujeme při každém načtení úvodní stránky zobrazující banner.

Vzhledem k tomu, že na začátku testování obvykle o CTR jednotlivých bannerů nevíme nic, algoritmus je u všech bannerů inicializován hodnotami α = β = 1. Žádný z nich tedy není preferovaný a všechny mají stejnou šanci na zobrazení. Dalšími kroky se odhad zpřesňuje a algoritmus postupně konverguje do stavu, kdy bude téměř vždy vyhrávat banner s nejvyšší skutečnou hodnotou CTR.

Jeden z okamžiků rozhodování mezi bannery k zobrazení je zachycen na obrázku 2.

Obrázek 2: Bayesovský bandita.

Obrázek 2: Bayesovský bandita.

Jedná se o stav, kdy už byl každý z bannerů několikrát zobrazen a máme tedy hrubý odhad CTR. Střední hodnota banneru 1 na obrázku je 0.25, u banneru 2 je 0.4 a u banneru 3 je 0.67. Nejvyšší pravděpodobnost zobrazení má tedy banner 3. Nemusí však nutně zvítězit, protože aktuální CTR losujeme z příslušných beta rozdělení a například oblast okolo 0.5 má nenulovou pravděpodobnost u všech tří bannerů. Svislé čárkované čáry představují příklady vylosovaných hodnot. Vidíme, že oproti očekávání byl banner 2 předstižen bannerem 1, avšak nejvyšší hodnoty dosáhl přece jen banner 3. Proto ho zobrazíme.

Jednoduchou implementaci algoritmu bayesovského bandity představuje následující kód v jazyce Python.


from scipy.stats import beta
from random import random

class Banner():
    def __init__(self, CTR, alpha=1, beta=1):
        self.CTR = CTR
        self.alpha = alpha
        self.beta = beta

    def update(self, click):
        if click == True:
            self.alpha += 1
        else:
            self.beta += 1

    def sample(self):
        return beta(self.alpha, self.beta).rvs()

    def getCTR(self):
        return self.CTR

class BayesBandit:
    def __init__(self, CTRs):
        self.banners = []
        for CTR in CTRs:
            self.banners.append(Banner(CTR))

    def selectBanner(self):
        sampledCTRs = [banner.sample() for banner in self.banners]
        return sampledCTRs.index(max(sampledCTRs))

    def simulateUser(self, banner):
        CTR = self.banners[banner].getCTR()
        if random() < CTR:
            self.banners[banner].update(False)
        else:
            self.banners[banner].update(True)

if __name__ == "__main__":
    bandit = BayesBandit([0.25, 0.4, 0.67])
    for i in range(100):
        banner = bandit.selectBanner()
        bandit.simulateUser(banner)
        print "Banner %d" % (banner + 1)

Třída Banner, definovaná na řádcích 4-20, implementuje jednotlivé bannery. Je inicializovaná skutečnou hodnotou CTR a volitelně počátečními parametry α a β. Metoda update provádí aktualizaci parametrů na základě toho, jestli uživatel klikl na banner poté, co se mu zobrazil. Metoda sample vybírá náhodnou hodnotu CTR z odpovídajícího beta rozdělení. Metoda selectBanner třídy BayesBandit nejprve vylosuje CTR pro všechny bannery a poté vrátí index banneru s nejvyšší hodnotou CTR. Simulace chování skutečného uživatele je prováděna v metodě simulateUser.  Ta využívá skutečné hodnoty CTR banneru a rozhoduje o tom, zda by skutečný uživatel na banner klikl nebo ne. Samotný testovací program je zapsaný na řádcích 39-44. Nejprve jsou vytvořeny 3 bannery se skutečnými hodnotami CTR 0.25, 0.4 a 0.67. Poté je opakovaně vybírán banner k zobrazení a simulováno chování uživatele. Číslo zobrazeného banneru je v každém kroku vypsáno na standardní výstup. Po spuštění programu je vidět, že algoritmus zpočátku vybírá různé bannery, ale s rostoucím počtem iterací začíná preferovat banner 3. Už po dosažení 100 iterací je zobrazován téměř výhradně banner s nejvyšším skutečným CTR.

Popsaný algoritmus předpokládá, že o skutečných hodnotách CTR před zahájením testování nic nevíme. To však nemusí platit vždy. Provozovatel e-shopu například může předem usoudit, že banner 3 bude u uživatelů oblíbenější než banner 1. Potom je možné nastavit vhodné iniciální hodnoty α a β namísto defaultních 1. Pokud bude odhad správný, konvergenci to ještě urychlí. Pokud by však byl iniciální expertní odhad výrazně odlišný od skutečných hodnot, algoritmus sice po čase dokonverguje do správných hodnot, ale bude k tomu zapotřebí mnohem více iterací.

Pravděpodobnostní programovací jazyky

Bayes

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í

N (x | \mu, \sigma^2) = \dfrac{1}{\sigma\sqrt{2\pi}} e^{\dfrac{-(x-\mu)^2}{2\sigma^2}}

jak je ilustrováno na obrázku 1.

Pravděpodobnostní rozdělení výšek lidí v České republice.

Obr. 1: Pravděpodobnostní rozdělení výšek lidí.

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ž, σ2)
  4. v případě, že jde o ženu, vyber její výšku z normálního rozdělení N(μžena, σ2)

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 &lt;- 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 &lt;- rnorm(5000, 180, 10)
zeny &lt;- rnorm(5000, 167, 10)
x &lt;- sample(c(muzi, zeny))
N &lt;- 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] &lt;- mu[1] + abs(theta)
    P[1] &lt;- 0.5
    P[2] &lt;- 0.5
    sigma ~ dunif(0, 50)
    tau &lt;- 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ě

y(x) = w_3x^3 + w_2x^2 +w_1x + w_0

a hledanými parametry koeficienty w = (w3, w2, w1, w0), jak je ilustrováno na obrázku 2.

 

Obr. 2: Regresní analýza.

Obr. 2: Regresní analýza.

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 &lt;- seq(-10, 10, by = 0.01)
N &lt;- length(x)
w3 &lt;- 1
w2 &lt;- 0
w1 &lt;- -1
w0 &lt;- 10
epsilon &lt;- rnorm(N, 0, 5)
y &lt;- 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í

y(x) = w_3x^3 - w_1x + 10

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] &lt;- 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 &lt;- 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 &lt;- 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