Laseme põhidefinitsiooni kõrvale jätta: Markovi keti Monte Carlo (MCMC) meetodid võimaldavad meil arvutada jaotuse proovid, kuigi me ei saa seda arvutada.

Mida see tähendab? Varundame ja räägime Monte Carlo proovivõtust.

Monte Carlo proovide võtmine

Mis on Monte Carlo meetodid?

'Monte Carlo meetodid või Monte Carlo katsed on lai arvutusalgoritmide klass, mis tugineb arvuliste tulemuste saamiseks korduvale juhuslikule valimisele.' ( Vikipeedia )

Lahustame selle.

Näide 1: Pindala hinnang

Kujutage ette, et teil on ebakorrapärane kuju, nagu allpool toodud kuju:

Ja teie ülesandeks on määrata selle kujuga ümbritsetud ala. Üks meetoditest, mida saate kasutada, on kujundisse väikeste ruutude tegemine, ruutude loendamine ja see annab teile ala üsna täpse ligikaudse ülevaate. See on aga keeruline ja aeganõudev.

Monte Carlo proovide võtmine päästmiseks!

Kõigepealt joonistame kuju ümber teadaoleva ala suure ruudu, näiteks 50 cm2. Nüüd 'riputame' selle ruudu ja hakkame viskama noolt suvaliselt kuju järgi.

Järgmisena loendame ristkülikus olevate noolemängude koguarvu ja meid huvitava kuju noolte arvu. Oletame, et kasutatud noolemänge oli kokku 100 ja 22 neist jõudsid kuju sisse. Nüüd saab pindala arvutada lihtsa valemi abil:

kuju pind = ruudu pindala * (noolel olevate noolte arv) / (noolte arv ruudul)

Nii et meie puhul taandub see:

kuju pind = 50 * 22/100 = 11 cm2

Kui korrutada “noolemängu” arv 10-kordse teguriga, saab see lähendus tegelikule vastusele väga lähedaseks:

kuju pind = 50 * 280/1000 = 14 cm2

Nii jagame keerulised ülesanded, nagu eespool toodud, Monte Carlo valimi abil.

Suurte arvude seadus

Ala lähendamine oli lähemal, seda rohkem nooli visati, ja see on tingitud Suurte arvude seadus :

“Suurte arvude seadus on lause, mis kirjeldab sama katse palju kordi sooritamise tulemust. Seaduse järgi peaks suure hulga katsete tulemuste keskmine olema oodatud väärtuse lähedal ja kipub lähemale jääma, kui rohkem katseid tehakse. ”

See viib meid järgmise näite, kuulsa Monty Halli probleemi juurde.

Näide 2: Monty Halli probleem

The Monty Halli probleem on väga kuulus ajutrust:

“Seal on kolm ust, ühel on auto taga, teistel kits taga. Valite ukse, peremees avab teise ukse ja näitab, et selle taga on kits. Seejärel küsib ta teilt, kas soovite oma otsust muuta. Kas sa? Miks? Miks mitte?'

Esimene asi, mis pähe tuleb, on see, et võiduvõimalused on võrdsed, olenemata sellest, kas vahetate või mitte, kuid see pole tõsi. Teeme sama demonstreerimiseks lihtsa vooskeemi.

Eeldades, et auto on ukse 3 taga:

Seega, kui vahetate, võidate ⅔ korda ja kui ei vahetata, siis ainult ⅓ korda.

Lahendame selle valimi abil.

wins = [] for i in range(int(10e6)): car_door = assign_car_door() choice = random.randint(0, 2) opened_door = assign_door_to_open(car_door) did_he_win = win_or_lose(choice, car_door, opened_door, switch = False) wins.append(did_he_win) print(sum(wins)/len(wins))

assign_car_door() funktsioon on lihtsalt juhuslike arvude generaator, mis valib ukse 0, 1 või 2, mille taga on auto. Kasutades assign_door_to_open valib ukse, mille taga on kits ja mis pole teie valitud, ning peremees avab selle. win_or_lose naaseb tõsi või vale , mis tähistab, kas olete auto võitnud või mitte, on vaja bool-lülitit, mis ütleb, kas vahetasite ukse või mitte.

Käivitame selle simulatsiooni 10 miljonit korda:

• Võidu tõenäosus, kui te ei vaheta: 0.334134
• Võidu tõenäosus, kui vahetate: 0.667255

See on väga lähedal vastustele, mida vooskeem meile andis.

Tegelikult, mida rohkem me seda simulatsiooni käivitame, seda lähemale jõuab vastus tegelikule väärtusele, kinnitades seega suurte arvude seadust:

Sama võib näha ka sellest tabelist:

Simulatsioonid töötavad	Võidu tõenäosus, kui vahetate	Võidu tõenäosus, kui te ei vaheta
10	0.6	0.2
10 ^ 2	0,63	0,33
10 ^ 3	0.661	0,333
10 ^ 4	0,6683	0,3236
10 ^ 5	0.66762	0,33171
10 ^ 6	0.667255	0,334134
10 ^ 7	0,6668756	0,3332821

Bayesi mõtteviis

'Frequentist, tuntud kui statistika klassikalisem versioon, eeldab, et tõenäosus on sündmuste pikaajaline sagedus (seega ka antud tiitel).'

„Bayesi statistika on statistikavaldkonna teooria, mis põhineb tõenäosuse Bayesi tõlgendusel, kus tõenäosus väljendab teatud usku mingisse sündmusesse. Uskumuse määr võib põhineda varasematel teadmistel sündmuse kohta, näiteks varasemate katsete tulemustel, või isiklikel veendumustel sündmuse kohta. ' - alates Tõenäosuslik programmeerimine ja häkkerite Bayesi meetodid

Mida see tähendab?

Fraktsistlikus mõtteviisis vaatleme tõenäosusi pikas perspektiivis. Kui sagedane esindaja ütleb, et autoõnnetuse tõenäosus on 0,001%, tähendab see, et kui arvestada lõpmatuid autoreise, lõpeb neist 0,001% krahhiga.

Bayesi mõtteviis on erinev, kuna alustame priori, veendumusest. Kui me räägime veendumusest 0, tähendab see, et teie veendumus on, et sündmus ei juhtu kunagi; vastupidi, uskumus 1 tähendab, et olete kindel, et see juhtub.

Seejärel, kui hakkame andmeid jälgima, ajakohastame seda veendumust, et neid arvesse võtta. Kuidas me seda teeme? Kasutades Bayesi teoreem .

Bayesi teoreem

.... (üks)

Jagame selle laiali.

• P(A | B) annab meile sündmuse A tõenäolise sündmuse tõenäosuse B. See on tagumine, B on andmed, mida me vaatasime, nii et me sisuliselt ütleme, milline on sündmuse toimumise tõenäosus, arvestades meie täheldatud andmeid.
• P(A) on prior, meie usk, et sündmus A juhtub.
• P(B | A) on tõenäosus, kui suur on tõenäosus, et me jälgime andmeid, arvestades, et A on tõene.

Vaatame näiteks vähi sõeluuringut.

Vähi tõenäosus

Oletame, et patsient läheb mammograafiat tegema ja mammograafia on positiivne. Kui suur on tõenäosus, et patsiendil on tegelikult vähk?

Määratleme tõenäosused:

• 1% naistest põeb rinnavähki.
• Mammogrammid tuvastavad vähi 80% juhtudest, kui see tegelikult esineb.
• 9,6% mammograafidest väidab valesti, et teil on vähk, kui teil seda tegelikult pole.

Seega, kui te ütleksite, et kui mammogramm peaks positiivse tähenduse saama, on 80% tõenäosus, et teil on vähk, oleks see vale. Te ei võtaks arvesse, et vähktõbi on haruldane sündmus, st ainult 1% naistest on rinnavähk. Me peame seda võtma varasemana ja siin tuleb mängu Bayesi teoreem:

P (C + | T +) = (P (T + | C +) * P (C +)) / P (T +)

• P(C+ | T+): See on tõenäosus, et vähk on olemas, arvestades, et test oli positiivne, see on see, mis meid huvitab.
• P(T+ | C+): See on tõenäosus, et test on positiivne, arvestades vähki, on see, nagu eespool arutletud, võrdne 80% = 0,8.
• P(C+): See on varasem tõenäosus, vähi põdeva inimese tõenäosus, mis võrdub 1% = 0,01.
• P(T+): See on tõenäosus, et test on positiivne, olenemata sellest, seega on sellel kaks komponenti: P (T +) = P (T + | C-) P (C -) + P (T + | C +) P (C +)
• P(T+ | C-): See on tõenäosus, et testi tulemus oli positiivne, kuid vähki pole, selle annab 0,096.
• P(C-): See on vähktõve tõenäosus, kuna vähi tõenäosus on 1%, see on võrdne 99% = 0,99.
• P(T+ | C+): See on tõenäosus, et testi tulemus oli positiivne, arvestades, et teil on vähk, see on 80% = 0,8.
• P(C+): See on vähktõve tõenäosus, mis võrdub 1% = 0,01.

Selle kõige ühendamine algsesse valemisse:

Seega, arvestades, et mammograafia tulemus oli positiivne, on 7,76% tõenäosus, et patsiendil on vähk. Esialgu võib see tunduda kummaline, kuid on mõistlik. Test annab valepositiivse 9,6% ajast (üsna kõrge), seega on antud populatsioonis palju valepositiivseid. Haruldase haiguse puhul on enamus positiivseid testitulemusi valed.

Vaadakem nüüd uuesti üle Monty Halli probleem ja lahendage see Bayesi teoreemi abil.

Monty Halli probleem

Priorid võib määratleda järgmiselt:

• Oletame, et valisite alguses ukse A.
• P (H) = ⅓, ⅓. ⅓ kõigi kolme ukse puhul, mis tähendab, et enne ukse avamist ja kitse paljastamist oli auto võrdne võimalus olla nende kõigi taga.

Tõenäosust saab määratleda järgmiselt:

• P(D|H), kus sündmus D on see, et Monty valib ukse B ja selle taga pole ühtegi autot.
• P(D|H) = ½, kui auto on ukse A taga - kuna on 50% tõenäosus, et ta valib ukse B, ja 50% tõenäosus, et ta valib ukse C
• P(D|H) = 0, kui auto on ukse B taga, kuna on 0% tõenäosus, et ta valib ukse B, kui auto on selle taga.
• P(D|H) = 1, kui auto on C taga ja valite A, on 100% tõenäosus, et ta valib ukse B

Meid huvitab P(H|D), mis on tõenäosus, et auto on ukse taga, arvestades, et see näitab meile kitse ühe ukse taga.

Tagantpoolt on P(H|D) näha, et A-le uksele C üleminek suurendab võidu tõenäosust.

Metropolis-Hastings

Ühendagem nüüd kõik, mida siiani käsitlesime, ja proovime mõista, kuidas Metropolis-Hastings töötab.

Metropolis-Hastings kasutab Bayesi teoreemi kompleksse jaotuse tagumise jaotuse saamiseks, millest valimi võtmine on keeruline.

Kuidas? Põhimõtteliselt valib see ruumist juhuslikult erinevad proovid ja kontrollib, kas uus proov on tõenäolisem tagumisest kui viimane proov, kuna vaatleme tõenäosuste suhet, P (B) võrrandis (1) saab tühistatud:

P (aktsepteerimine) = (P (uus proov) * uue valimi tõenäosus) / (P (vana proov) * vana valimi tõenäosus)

Iga uue valimi tõenäosuse otsustab funktsioon f . Sellepärast f peab olema proportsionaalne tagumisega, millest proovime proovida.

Otsustamaks, kas θ ′ aktsepteeritakse või lükatakse tagasi, tuleb iga uue kavandatud θ jaoks arvutada järgmine suhe:

Kus θ on vana proov, P(D| θ) on valimi tõenäosus θ.

Kasutame näite, et sellest paremini aru saada. Oletame, et teil on andmeid, kuid soovite teada saada neile sobiva normaaljaotuse, nii et:

X ~ N (keskmine, standard)

Nüüd peame määratlema proorid keskmise ja std jaoks. Lihtsuse huvides eeldame, et mõlemad järgivad normaaljaotust keskmise 1 ja standardiga 2:

Keskmine ~ N (1, 2)
standard ~ N (1, 2)

Nüüd genereerime mõned andmed ja oletame, et tegelik keskmine ja standardväärtus on vastavalt 0,5 ja 1,2.

import numpy as np import matplotlib.pyplot as plt import scipy.stats as st meanX = 1.5 stdX = 1.2 X = np.random.normal(meanX, stdX, size = 1000) _ = plt.hist(X, bins = 50)

PyMC3

Kasutagem kõigepealt teeki nimega PyMC3 andmete modelleerimiseks ja keskmise ja standardväärtuse tagumise jaotuse leidmiseks.

import pymc3 as pm with pm.Model() as model: mean = pm.Normal('mean', mu = 1, sigma = 2) std = pm.Normal('std', mu = 1, sigma = 2) x = pm.Normal('X', mu = mean, sigma = std, observed = X) step = pm.Metropolis() trace = pm.sample(5000, step = step)

Läheme koodist üle.

Esiteks määratleme priori keskmise ja std jaoks, st norm, mille keskmine on 1 ja std 2.

x = pm. Normaalne (“X”, mu = keskmine, sigma = std, vaadeldud = X)

Selles reas määratleme meid huvitava muutuja; selleks on vaja varem määratletud keskmist ja standardväärtust, samuti tuleb arvestada varem arvutatud väärtusi.

Vaatame tulemusi:

_ = plt.hist(trace['std'], bins = 50, label = 'std') _ = plt.hist(trace['mean'], bins = 50, label = 'mean') _ = plt.legend()

Standard on keskmiselt umbes 1,2 ja keskmine on 1,55 - üsna lähedal!

Teeme seda nüüd nullist, et Metropolis-Hastingsist aru saada.

Algusest

Ettepaneku genereerimine

Kõigepealt mõistame, mida Metropolis-Hastings teeb. Selle artikli varem ütlesime, 'Metropolis-Hastings valib kosmosest juhuslikult erinevad proovid,' aga kuidas ta teab, millise valimi valida?

Ta teeb seda ettepanekute jaotuse abil. See on normaalne jaotus, mis on koondatud praegu aktsepteeritud valimisse ja mille STD on 0,5. Kui 0,5 on hüperparameeter, võime seda näpistada; mida rohkem on ettepaneku STD, seda kaugemale see praegu aktsepteeritud valimist otsib. Kodeerime selle.

Kuna peame leidma jaotuse keskmise ja standardi, võtab see funktsioon praegu aktsepteeritud keskmise ja praegu aktsepteeritud keskmise ning tagastab ettepanekud mõlema jaoks.

def get_proposal(mean_current, std_current, proposal_width = 0.5): return np.random.normal(mean_current, proposal_width), np.random.normal(std_current, proposal_width)

Ettepaneku vastuvõtmine / tagasilükkamine

Kodeerime nüüd loogika, mis ettepaneku aktsepteerib või tagasi lükkab. Sellel on kaks osa: prior ja tõenäosus .

Kõigepealt arvutame välja tõenäosuse, et ettepanek tuleb varasemast:

def prior(mean, std, prior_mean, prior_std): return st.norm(prior_mean[0], prior_mean[1]).pdf(mean)* st.norm(prior_std[0], prior_std[1]).pdf(std)

See lihtsalt võtab tähendab ja tundi ja arvutab, kui tõenäoline see on tähendab ja tundi tuli eelnev levitamine mida me määratlesime.

Tõenäosuse arvutamisel arvutame välja, kui tõenäoline on, et andmed, mida nägime, tulid kavandatud jaotusest.

def likelihood(mean, std, data): return np.prod(st.norm(mean, std).pdf(X))

Niisiis korrutame iga andmepunkti jaoks selle andmepunkti tõenäosuse kavandatud jaotusest.

Nüüd peame neid funktsioone kutsuma praeguse keskmise ja standardväärtuse ning pakutava jaoks tähendab ja tundi .

prior_current = prior(mean_current, std_current, prior_mean, prior_std) likelihood_current = likelihood(mean_current, std_current, data) prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std) likelihood_proposal = likelihood(mean_proposal, std_proposal, data)

Kogu kood:

def accept_proposal(mean_proposal, std_proposal, mean_current, std_current, prior_mean, prior_std, data): def prior(mean, std, prior_mean, prior_std): return st.norm(prior_mean[0], prior_mean[1]).pdf(mean)* st.norm(prior_std[0], prior_std[1]).pdf(std) def likelihood(mean, std, data): return np.prod(st.norm(mean, std).pdf(X)) prior_current = prior(mean_current, std_current, prior_mean, prior_std) likelihood_current = likelihood(mean_current, std_current, data) prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std) likelihood_proposal = likelihood(mean_proposal, std_proposal, data) return (prior_proposal * likelihood_proposal) / (prior_current * likelihood_current)

Nüüd loome lõpliku funktsiooni, mis seob kõik omavahel ja genereerib tagumised proovid, mida vajame.

Tagumise loomine

Nüüd nimetame ülalkirjutatud funktsioone ja genereerime tagumise!

def get_trace(data, samples = 5000): mean_prior = 1 std_prior = 2 mean_current = mean_prior std_current = std_prior trace = { 'mean': [mean_current], 'std': [std_current] } for i in tqdm(range(samples)): mean_proposal, std_proposal = get_proposal(mean_current, std_current) acceptance_prob = accept_proposal(mean_proposal, std_proposal, mean_current, std_current, [mean_prior, std_prior], [mean_prior, std_prior], data) if np.random.rand()

Arenguruum

Log on teie sõber!

Meenuta a * b = antilog(log(a) + log(b)) ja a / b = antilog(log(a) - log(b)).

See on meile kasulik, kuna tegeleme väga väikeste tõenäosustega, mille korrutamisel saadakse null, seega lisame pigem logiprobleemi, probleem lahendatud!

Niisiis saab meie eelmine kood:

def accept_proposal(mean_proposal, std_proposal, mean_current, std_current, prior_mean, prior_std, data): def prior(mean, std, prior_mean, prior_std): return st.norm(prior_mean[0], prior_mean[1]).logpdf(mean)+ st.norm(prior_std[0], prior_std[1]).logpdf(std) def likelihood(mean, std, data): return np.sum(st.norm(mean, std).logpdf(X)) prior_current = prior(mean_current, std_current, prior_mean, prior_std) likelihood_current = likelihood(mean_current, std_current, data) prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std) likelihood_proposal = likelihood(mean_proposal, std_proposal, data) return (prior_proposal + likelihood_proposal) - (prior_current + likelihood_current)

Kuna tagastame aktsepteerimise tõenäosuse logi:

if np.random.rand()

Saab:

if math.log(np.random.rand())

Käivitame uue koodi ja koostame tulemused.

_ = plt.hist(trace['std'], bins = 50, label = 'std') _ = plt.hist(trace['mean'], bins = 50, label = 'mean') _ = plt.legend()

Nagu näete, tundi on keskendatud punktile 1.2 ja tähendab kell 1.5. Me tegime seda!

Edasi liikuma

Nüüdseks saate loodetavasti aru, kuidas Metropolis-Hastings töötab, ja võite mõelda, kus saaksite seda kasutada.

Noh, Bayesi meetodid häkkeritele on suurepärane raamat, mis selgitab tõenäosuslikku programmeerimist ja kui soovite rohkem teada saada Bayesi teoreemi ja selle rakenduste kohta, Mõelge Bayesile on Allen B. Downey suurepärane raamat.

Täname lugemast ja loodan, et see artikkel julgustab teid avastama Bayesi statistika imelist maailma.

Põhitõdesid mõistes

Mida tähistab MCMC?

MCMC tähistab Markovi ketti Monte Carlo, mis on tõenäosusjaotise valimi moodustamise meetodite klass.

Milline on Bayesi lähenemine otsuste tegemisele?

Bayesi lähenemisviisis otsuste tegemisel alustate kõigepealt prioriteedist, see on teie veendumused, siis kui andmed tulevad, lisate need andmed nende priorite värskendamiseks tagumise poole saamiseks.

Mis on Bayesi mudel?

Bayesi mudel on statistiline mudel, kus kogu mudeli ebakindluse esitamiseks kasutate tõenäosust.