Lanțuri Markov Monte Carlo (MCMC): Eșantionarea distribuțiilor de probabilitate complexe

Unele distribuții de probabilitate sunt pur și simplu prea complexe ca să poți lucra direct cu ele.

Când modelezi date din lumea reală, matematica cedează adesea înainte să ajungi la ceva util. De multe ori, integralele par gestionabile pe hârtie, dar devin intratabile în momentul în care adaugi câteva variabile latente. Asta e mai ales frecvent în inferența bayesiană, unde distribuția a posteriori combină convingerile anterioare cu datele observate – iar rezultatul e ceva ce nu poți rezuma printr-o formulă simplă.

Ideea de bază din Markov Chain Monte Carlo este că, în loc să ataci direct matematica, MCMC explorează distribuția prin simulare. Extrage eșantioane care îi reflectă forma fără să fie nevoie s-o calculezi complet.

În acest articol, voi acoperi conceptele de bază din spatele MCMC, voi parcurge cele mai comune algoritmi și îți voi arăta cum să îl pui la treabă în Python.

Ai nevoie de un refresh la matematică în Python? Citește articolele noastre de blog Demystifying Mathematical Concepts for Deep Learning ca să vezi matematica aplicată în Numpy.

Ce este Markov Chain Monte Carlo?

Markov Chain Monte Carlo (MCMC) este o familie de algoritmi care generează eșantioane din distribuții de probabilitate – chiar și atunci când acele distribuții sunt prea complexe pentru a lucra direct cu ele.

Numele are două părți. Lanțul Markov controlează cum se mișcă algoritmul prin stările posibile. Fiecare pas depinde doar de unde ești acum, nu de întreaga istorie a modului în care ai ajuns acolo. Partea Monte Carlo înseamnă că folosești eșantionare aleatorie pentru a estima cantitățile de interes.

Împreună, MCMC construiește un lanț de eșantioane aleatorii care, în timp, reflectă forma distribuției tale țintă. Este, în primul rând, o tehnică de eșantionare. Nu rezolvi matematica exact, ci o aproximezi prin simulare.

De ce este nevoie de Markov Chain Monte Carlo

Problema cu distribuțiile datelor din lumea reală este că nu sunt nici pe departe la fel de curate ca cele din manuale.

În inferența bayesiană, încerci adesea să calculezi o distribuție a posteriori – probabilitatea actualizată a parametrilor modelului după ce ai văzut datele. Formula pare simplă pe hârtie: înmulțești priorul cu verosimilitatea, apoi împarți la verosimilitatea marginală. Ultimul termen cere integrarea peste fiecare valoare posibilă a parametrilor. În dimensiuni mari, acea integrală este practic imposibil de calculat.

Devine și mai rău pe măsură ce modelul crește. Când adaugi mai mulți parametri sau mai multe variabile latente, calculul exact devine un drum înfundat. Te vei lovi de asta într-o serie de situații comune:

• Modele ierarhice bayesiene: Unde parametrii depind de alți parametri, creând distribuții imbricate care nu se integrează curat
• Posteriori de înaltă dimensiune: Unde numărul de parametri face ca aproximațiile pe grilă să fie imposibil de realizat computațional
• Prioruri neconjugate: Unde priorul și verosimilitatea nu se combină într-o distribuție recognoscibilă

MCMC este o soluție bună în aceste scenarii. În loc să calculeze distribuția, extrage eșantioane din ea. Aceste eșantioane conțin tot ce ai nevoie fără să rezolvi integralul.

Cele două idei din spatele MCMC

MCMC combină două idei care, luate separat, sunt destul de simple, dar puternice împreună. Hai să ți le explic.

Lanțuri Markov

Un lanț Markov este o succesiune de stări în care fiecare pas depinde doar de locul în care te afli acum.

Unde ai fost înainte nu contează. Doar starea curentă determină unde mergi mai departe. Această proprietate „fără memorie” – numită formal proprietatea Markov – este ceea ce face matematica ușor de gestionat și algoritmul practic.

Lanțul se deplasează prin stări pas cu pas și, cu setarea potrivită, în cele din urmă se stabilizează într-o distribuție staționară – un tipar stabil în care probabilitatea de a fi într-o anumită stare nu mai se schimbă. Acea distribuție stabilă este exact ceea ce MCMC este conceput să folosească.

Metode Monte Carlo

Metodele Monte Carlo folosesc eșantionarea aleatorie pentru a estima cantități greu de calculat direct.

Ideea este să extragi suficiente eșantioane aleatorii dintr-o distribuție, apoi să îi estimezi media, varianța sau orice altă proprietate doar uitându-te la eșantioane. Cu cât extragi mai multe eșantioane, cu atât estimările se apropie mai mult de valorile reale.

Luate separat, metodele Monte Carlo cer să eșantionezi direct din distribuție – ceea ce este exact problema pe care încercăm s-o rezolvăm. Lanțurile Markov se ocupă de acea parte.

Cum funcționează Markov Chain Monte Carlo

MCMC este o buclă cu o decizie simplă la fiecare pas.

1. Pornește de la o stare inițială: Alege un punct de start în spațiul parametrilor, chiar dacă este o aproximare grosieră
2. Propune o stare nouă: Folosește un mecanism de propunere pentru a sugera unde să te muți în continuare
3. Decide dacă o accepți: Compară starea propusă cu cea curentă și acceptă sau respinge în funcție de cât de probabilă este fiecare sub distribuția țintă
4. Repetă: Rulează această buclă de mii de ori, colectând eșantioane pe parcurs
5. Folosește eșantioanele: Estimează medii, intervale credibile sau orice altă cantitate de interes din eșantioanele colectate

Pasul de acceptare/respingere este locul unde se întâmplă „magia”.

Acceptând stări mai bune mai des decât pe cele mai slabe, lanțul gravitează către regiuni cu probabilitate mare – fără să fie nevoie să calculeze distribuția completă.

Eșantioanele timpurii depind de locul din care ai pornit, așa că acelea sunt aruncate. După suficiente iterații, lanțul își „uită” punctul de plecare, iar eșantioanele rămase reflectă forma reală a distribuției țintă.

Distribuția țintă și distribuția staționară

MCMC este construit în jurul obiectivului de a genera eșantioane dintr-o distribuție țintă din care nu poți eșantiona direct.

Distribuția țintă este ceea ce încerci să înțelegi – de obicei o distribuție a posteriori în inferența bayesiană. Îi cunoști forma până la o constantă de normalizare, dar nu poți calcula direct acea constantă. MCMC nu are nevoie de ea.

Fiecare algoritm MCMC este conceput astfel încât lanțul său Markov să aibă o distribuție staționară care să se potrivească cu distribuția țintă. O distribuție staționară este distribuția în care lanțul se stabilizează după suficiente pași.

Lasă lanțul să ruleze și va începe să producă eșantioane care arată exact ca extrageri din distribuția ta țintă. Integralul este ocolit.

Algoritmi populari Markov Chain Monte Carlo

Există câțiva algoritmi MCMC pe care îi vei întâlni în practică. Toți urmează aceeași buclă de bază, dar diferă prin modul în care propun stări noi și prin modul în care folosesc informațiile despre distribuția țintă.

Algoritmul Metropolis

Algoritmul Metropolis este cel mai simplu algoritm MCMC și cel care a pornit totul.

La fiecare pas, propune o stare nouă adăugând zgomot aleator la cea curentă. Dacă starea propusă are o probabilitate mai mare sub distribuția țintă, este acceptată întotdeauna. Dacă e mai mică, este acceptată cu o probabilitate proporțională cu raportul dintre cele două probabilități – altfel, lanțul rămâne pe loc.

Acest mecanism de acceptare/respingere înseamnă că lanțul va petrece mai mult timp în regiunile cu probabilitate mare fără să calculeze vreodată distribuția completă.

Algoritmul Metropolis folosește o distribuție de propunere simetrică, ceea ce înseamnă că are șanse egale să propună un pas în orice direcție. Tinde să cedeze pe măsură ce modelele cresc.

Algoritmul Metropolis-Hastings

Algoritmul Metropolis-Hastings (MH) generalizează Metropolis permițând distribuții de propunere asimetrice.

MH ajustează probabilitatea de acceptare pentru a ține cont de faptul că unele propuneri sunt mai probabile decât altele. Poți regla propunerea în funcție de forma distribuției țintă, ceea ce duce la o explorare mai bună și la o convergență mai rapidă.

Majoritatea metodelor MCMC moderne sunt extensii ale MH sau construite pe aceleași principii. Așadar, dacă înțelegi Metropolis-Hastings, înțelegi fundația domeniului.

Eșantionarea Gibbs

Eșantionarea Gibbs actualizează câte o variabilă pe rând, în loc să propună o stare nouă pentru toți parametrii simultan.

La fiecare pas, eșantionează fiecare variabilă din distribuția sa condițională – distribuția acelei variabile date valorile curente ale tuturor celorlalte. Când ai parcurs toate variabilele, ai completat o iterație întreagă.

Aceasta evită complet pasul de acceptare/respingere, deoarece fiecare extragere condițională este acceptată întotdeauna. E utilă când distribuția comună completă e greu de eșantionat, dar condiționalele sunt tratabile, lucru des întâlnit în modelele ierarhice bayesiene.

Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) este primul algoritm care a făcut inferența bayesiană modernă practică la scară.

În loc să propună aleatoriu stări noi, HMC folosește informația de gradient din distribuția țintă pentru a propune stări aflate departe de poziția curentă, dar care au totuși șanse mari să fie acceptate. Se mișcă mult mai bine prin spațiul parametrilor decât metodele de tip random-walk. Există mai puține propuneri respinse și o explorare mai bună a distribuțiilor de înaltă dimensiune.

Metodele random-walk precum Metropolis nu scalează pe măsură ce crește numărul parametrilor. HMC nu are această problemă în aceeași măsură.

HMC este motorul din spatele Stan, una dintre cele mai folosite platforme de programare probabilistică. No-U-Turn Sampler (NUTS), o extensie adaptivă a HMC folosită în PyMC, elimină nevoia de a regla manual dimensiunea pasului și numărul de pași.

Markov Chain Monte Carlo în statistica bayesiană

Dacă există un domeniu în care MCMC a avut cel mai mare impact, acela este inferența bayesiană.

Statistica bayesiană se concentrează pe distribuția a posteriori, adică probabilitatea actualizată a parametrilor modelului după observarea datelor. Calcularea ei înseamnă înmulțirea priorului cu verosimilitatea și normalizarea. Pasul de normalizare cere o integrală care rareori este tratabilă.

MCMC elimină complet acest pas. Evaluezi doar a posterioriul nenormalizat într-un punct dat și lași lanțul să facă restul.

Iată un exemplu simplu. Presupune că estimezi biasul unei monede. Pornești de la convingerea anterioară că moneda este probabil echitabilă, apoi observi o secvență de aruncări. Pentru un model simplu de monedă, a posterioriul are formă închisă. Dacă adaugi o structură ierarhică, adică estimezi biasul pentru o sută de monede simultan, devine imposibil de calculat.

Cu MCMC, rulezi lanțul, colectezi eșantioane din a posteriori și folosești acele eșantioane pentru a calcula ce ai nevoie.

Înțelegerea burn-in, convergenței și amestecării (mixing)

Aceste trei concepte îi încurcă pe data scientist-ii noi în MCMC. Dacă le înțelegi greșit, vei obține rezultate, dar nu vei ști de ce nu sunt de încredere.

Burn-in

Când un lanț Markov pornește, nu are idee unde sunt regiunile cu probabilitate mare ale distribuției țintă.

Eșantioanele timpurii sunt influențate de punctul tău de pornire, nu de distribuția țintă. Burn-in este practica de a le arunca. Rulezi lanțul un anumit număr de iterații la început, arunci acele eșantioane și păstrezi doar ce vine după ce lanțul a avut timp să găsească un punct de pornire bun.

Nu există o regulă universală pentru cât ar trebui să fie burn-in-ul. Depinde de model, de punctul de pornire și de cât de bine se amestecă lanțul. În practică, îl diagnostichezi vizual cu grafice de trasă (trace plots), nu alegând un număr fix din start.

Convergență

Convergența înseamnă că lanțul nu mai este influențat de punctul de pornire și acum extrage eșantioane care reflectă distribuția țintă.

Un lanț care nu a convergat produce eșantioane părtinitoare. Media pe care o calculezi din ele nu va corespunde mediei reale a posteriori. În schimb, va reflecta zona în care lanțul a rămas blocat.

Convergența este ceva ce evaluezi ulterior, cu diagnostice. Rularea mai multor lanțuri din puncte de pornire diferite și verificarea dacă sunt de acord este una dintre cele mai fiabile metode de a depista eșecurile de convergență.

Amestecare (mixing)

Un lanț care converge dar se amestecă prost este tot problematic.

Amestecarea descrie cât de bine explorează lanțul distribuția țintă. Un lanț care se amestecă bine se mișcă liber, vizitând atât regiuni cu probabilitate mare, cât și mică, și produce eșantioane aproximativ independente unele de altele. Un lanț care se amestecă prost rămâne într-o regiune multe iterații înainte să se miște și produce eșantioane puternic corelate care nu reprezintă întreaga distribuție.

Amestecarea slabă apare adesea ca un grafic de trasă care arată ca un râu lent, șerpuitor, în loc de o bandă orizontală zgomotoasă. Când vezi asta, e nevoie să reglezi samplerul – o distribuție de propunere mai bună sau un algoritm complet diferit.

Grafic de comparație a amestecării

Cum să evaluezi rezultatele MCMC

Îți voi arăta acum patru moduri în care poți evalua MCMC și când să folosești fiecare metodă.

Grafice de trasă (trace plots)

Un grafic de trasă arată valoarea eșantionată a unui parametru la fiecare iterație. Este primul lucru pe care ar trebui să-l verifici după ce rulezi MCMC.

Un grafic de trasă „sănătos” arată ca un zgomot alb în jurul unei medii stabile. N-ar trebui să vezi trenduri, porțiuni lungi plate sau derive lente. Dacă vezi lanțul rătăcind sau blocându-se într-o regiune multe iterații, e o problemă de amestecare și eșantioanele tale nu sunt de încredere.

Grafice de trasă vizualizate

Autocorelație

Eșantioanele MCMC nu sunt niciodată complet independente. Fiecare eșantion este influențat de cel dinainte. Autocorelația măsoară cât de puternic sunt corelate eșantioanele de-a lungul iterațiilor.

Autocorelația mare înseamnă că eșantioanele tale poartă mai puțină informație decât sugerează numărul lor. Două mii de eșantioane corelate îți pot oferi aceeași informație ca două sute independente. Majoritatea bibliotecilor MCMC includ grafice de autocorelație, astfel încât să poți vedea cât de repede scade corelația pe măsură ce eșantioanele sunt mai depărtate.

Grafice de autocorelație vizualizate

Mărimea efectivă a eșantionului

Mărimea efectivă a eșantionului (ESS) traduce acea autocorelație într-un număr practic: câte eșantioane independente este echivalent lanțul tău.

Dacă ai extras 5.000 de eșantioane, dar ESS este 200, lucrezi cu puterea statistică a 200 de extrageri independente. Un ESS mic înseamnă că trebuie să rulezi lanțul mai mult, să reglezi samplerul sau ambele. Majoritatea practicienilor țintesc un ESS de cel puțin câteva sute per parametru înainte de a avea încredere în estimări.

Diagnostic de convergență

Când rulezi mai multe lanțuri, poți testa formal dacă au convergat la aceeași distribuție. Diagnosticul Gelman–Rubin, raportat ca R-hat, compară varianța din fiecare lanț cu varianța între lanțuri.

Un R-hat aproape de 1,0 înseamnă că lanțurile sunt de acord, ceea ce e un semn bun. Valori peste 1,01 sau 1,05 (în funcție de pragul folosit de biblioteca ta) sugerează că lanțurile n-au convergat și ai nevoie de mai multe iterații. Majoritatea bibliotecilor moderne, precum PyMC, calculează automat R-hat și afișează avertismente când este prea mare.

Markov Chain Monte Carlo în Python

Python are câteva biblioteci pentru MCMC, fiecare cu o filosofie diferită.

• PyMC: Cea mai abordabilă opțiune pentru majoritatea data scientist-ilor. Îți definești modelul în sintaxă Python, iar PyMC se ocupă de eșantionare în „spate” folosind NUTS (o versiune adaptivă a HMC). Este soluția de bază pentru modelare bayesiană în Python.
• CmdStanPy / PyStan: Interfețe Python către Stan, un limbaj dedicat de programare probabilistică. Stan compilează modelul tău în C++ înainte să-l ruleze, ceea ce îl face rapid și potrivit pentru modele ierarhice complexe. Compromisul este o curbă de învățare mai abruptă.
• NumPy: Fără bibliotecă MCMC dedicată, dar poți implementa de la zero algoritmi precum Metropolis-Hastings. Util pentru a înțelege ce se întâmplă „sub capotă” înainte de a trece la instrumente de nivel mai înalt.

Pentru munca practică, PyMC este locul de pornire. Asta voi folosi, așa că, dacă lucrezi în paralel, asigură-te că instalezi mai întâi biblioteca:

pip install pymc

Ca să păstrez lucrurile simple, voi rămâne la un exemplu ușor: estimarea biasului unei monede dintr-o secvență de aruncări.

Pasul 1: Definește modelul

import pymc as pm
import numpy as np

# 1 = heads, 0 = tails
observed_flips = np.array([1, 0, 1, 1, 0, 1, 1, 1, 0, 1])

with pm.Model() as coin_model:
    # Prior: we believe the coin is probably fair
    bias = pm.Beta("bias", alpha=2, beta=2)

    # Likelihood: observed flips given the bias
    flips = pm.Bernoulli("flips", p=bias, observed=observed_flips)

Priorul pm.Beta codifică o credință slabă că moneda este echitabilă. Verosimilitatea pm.Bernoulli conectează modelul la datele observate.

Pasul 2: Rulează samplerul

with coin_model:
    trace = pm.sample(2000, tune=1000, return_inferencedata=True)

Rezultatul rulării samplerului

tune setează numărul de pași de burn-in – acele eșantioane sunt aruncate. sample extrage 2000 de eșantioane a posteriori per lanț după etapă de tuning.

Pasul 3: Inspectează a posterioriul

import arviz as az

az.plot_trace(trace, var_names=["bias"])
az.summary(trace, var_names=["bias"])

Graficul de trasă al modelului și rezultatele sumarului

az.summary() îți oferă media a posteriori, deviația standard și R-hat pentru fiecare parametru. Dacă R-hat este aproape de 1,0, lanțurile au convergat. az.plot_trace() trasează graficul de trasă și distribuția a posteriori alăturat pentru fiecare parametru.

Pentru acest set de date – 7 capete din 10 aruncări – media a posteriori este 0,642 cu o deviație standard de 0,124. Aceasta reflectă evidența din date, rămânând aproape de priorul de monedă echitabilă. R-hat este 1,00 și ESS este mult peste 2000, deci lanțurile au convergat și eșantioanele sunt de încredere.

Greșeli comune cu MCMC

MCMC e ușor de rulat, dar la fel de ușor de folosit greșit. Iată greșelile care apar cel mai des.

• Presupunerea convergenței prea devreme: Doar pentru că samplerul a rulat fără erori nu înseamnă că a convergat. Verifică mereu R-hat și graficele de trasă înainte să ai încredere în rezultate. Un lanț poate părea stabil și totuși să fie blocat într-o regiune a spațiului parametrilor.
• Ignorarea burn-in-ului: Eșantioanele timpurii sunt părtinitoare către punctul tău de pornire, nu către distribuția țintă. Aruncă-le. Majoritatea bibliotecilor gestionează asta automat prin parametrul tune, dar verifică să nu incluzi din greșeală acele eșantioane în analiză.
• Amestecare slabă: Un lanț care se mișcă prea încet produce eșantioane puternic corelate, care poartă mai puțină informație decât sugerează numărul lor. Dacă graficul de trasă arată ca o derivă lentă, în loc de o bandă orizontală zgomotoasă, distribuția de propunere are nevoie de reglaj sau îți trebuie un algoritm mai bun. Trecerea de la Metropolis la NUTS rezolvă adesea problema direct.
• Prea puține eșantioane: Un ESS mic înseamnă că estimările tale nu sunt de încredere, chiar dacă lanțul a convergat tehnic. Dacă ESS este o fracție mică din totalul extragerilor, rulează lanțul mai mult sau îmbunătățește amestecarea înainte de a trage concluzii.

MCMC vs alte metode de eșantionare

MCMC nu este singura modalitate de a aproxima o distribuție. Iată cum se compară cu alternativele.

• Eșantionarea prin respingere extrage propuneri dintr-o distribuție mai simplă și le acceptă în funcție de cât de bine se potrivesc țintei. Este exactă atunci când funcționează, dar ratele de acceptare se prăbușesc în dimensiuni mari.
• Eșantionarea prin importanță ponderă eșantioanele dintr-o distribuție de propunere pentru a aproxima așteptările sub distribuția țintă. Este rapidă și funcționează bine în dimensiuni mici, dar varianța ponderilor explodează pe măsură ce crește dimensionalitatea. În practică, câteva eșantioane ajung să poarte aproape toată greutatea, ceea ce face estimările nesigure.
• Inferența variațională (VI) potrivește o distribuție mai simplă pentru a aproxima ținta, minimizând o măsură de divergență. VI este mult mai rapidă decât MCMC și scalează la seturi mari de date, dar introduce eroare de aproximație. Distribuția ajustată poate rata structura a posteriori pe care MCMC ar surprinde-o.

Pe scurt:

MCMC versus alternative

MCMC este alegerea potrivită când acuratețea contează mai mult decât viteza. Dacă ai nevoie să scalezi la seturi mari de date sau să rulezi inferență în timp real, inferența variațională ar putea merita compromisul de acuratețe.

Concluzie

MCMC este unul dintre acele instrumente care par intimidante din exterior, dar capătă sens odată ce înțelegi ce face de fapt – construiește un lanț de eșantioane care reflectă treptat forma unei distribuții pe care nu o poți calcula direct.

Este, de asemenea, mult mai ușor de înțeles dacă îl desfaci în părți: lanțuri Markov și metode Monte Carlo.

Rolul său în statistica bayesiană este greu de supraestimat. Distribuțiile a posteriori care altfel ar fi inaccesibile devin rezolvabile în momentul în care ai un sampler de încredere. De aceea MCMC este în centrul bibliotecilor de programare probabilistică precum PyMC și Stan.

Dar înainte de a trece la implementare, ar trebui să-ți fixezi intuiția. Înțelege de ce lanțul are nevoie de burn-in, ce înseamnă de fapt amestecarea și cum să citești un grafic de trasă. Codul în sine este partea ușoară, deoarece bibliotecile Python ascund toate abstracțiile în spatele unor apeluri de funcții simple.

Dacă vrei să devii priceput în machine learning, înscrie-te la Machine Learning Scientist in Python. 85 de ore de materiale te vor ajuta să fii pregătit pentru job în 2026.