(Acest articol a fost publicat pentru prima dată pe https://pacha.dev/blogși cu amabilitate a contribuit la R-bloggeri). (Puteți raporta problema legată de conținutul acestei pagini aici)
Doriți să vă distribuiți conținutul pe R-bloggeri? dați clic aici dacă aveți un blog, sau aici dacă nu aveți.
Înainte de conținutul principal: creez o comunitate R pe Grupuri Google. Vă puteți alătura grupului folosind acest formular.
Algoritmul Expectation-Maximization (EM) este un cadru de optimizare iterativ folosit pentru a găsi probabilitate maximă estimări ale parametrilor atunci când un model depinde de variabile neobservate, latente.
Pentru regresia liniară, algoritmul EM nu este necesar, deoarece Pătratele minime obișnuite (OLS) furnizează o soluție analitică în formă închisă, care este (hat{beta} = (X^TX)^{-1} (X^T y)). Cu toate acestea, încadrarea regresiei liniare prin algoritmul EM este un exercițiu destul de clarificator înainte de a trece la regresia Poisson (sau binom/logit), modele în care lipsesc datele (de exemplu, Tobit) sau în care setul de date este generat de un amestec de regresii liniare (linii multiple ascunse).
Pentru a înțelege cum se aplică EM, luați în considerare un set de date cu observații (n) și variabile (sau caracteristici) (p < n). Fiecare observație (i) corespunde unui vector (x_i = (x_{i1}, x_{i2}, ldots, x_{ip})^T in mathbb{R}^p), unui răspuns scalar (y_i), iar modelul este
( y_i = sum_{j = 1}^p beta_j x_{ij} + e_i,quad e_i sim N(0, sigma^2), )
unde (beta = (beta_1, beta_2, ldots, beta_p)^T in mathbb{R}^p) sunt ponderile care trebuie estimate.
Stivuirea tuturor observațiilor (n) dă forma matricei compacte (y = Xbeta + e), unde
( X = begin{pmatrix} x_{11} & x_{12} & ldots & x_{1p} \ x_{21} & x_{22} & ldots & x_{2p} \ vdots & ddots & & \ x_{n1} & x_{n2} & x_{n2} & ldots & l{ntrix} mathbb{R}^{n times p}, quad y = begin{pmatrix} y_1 \ y_2 \ vdots \ y_n end{pmatrix} in mathbb{R}^n, quad e = begin{pmatrix} e_1 \ e_2 \ e_2 \m e vdots {pmatrix} {0} sigma^2 I_n).
În cadrul EM tratăm variabilele neobservate (latente) (z_i = sum_{j = 1}^p beta_j x_{ij}) ca date complete, chiar dacă putem calcula distribuția lor exact. Ideea este de a stabili mașina care se generalizează la alte modele.
Configurare: probabilitatea de înregistrare a datelor complete
Dacă am observa atât (y_i) cât și (z_i), probabilitatea de înregistrare a datelor complete pentru (beta) și (sigma^2) ar fi
( ell_c(beta, sigma^2) = -frac{n}{2}log(2pisigma^2) – frac{1}{2sigma^2}sum_{i=1}^n (y_i – z_i)^2. )
deoarece (y_i mid z_i sim N(z_i, sigma^2)) (modelul de zgomot), iar (z_i) este determinist dat (beta).
Pas de așteptare
Având în vedere estimările actuale ale parametrilor (beta^{
Pentru regresia liniară variabila latentă (z_i) este complet determinată de (beta). Nu există nicio incertitudine odată ce (beta) este fixat, deci așteptarea condiționată este doar valoarea ajustată curentă
( mathbb{E}left(z_i mid y_i, beta^{
În notația matriceală, aceasta este (i)-a intrare a lui (X beta^{
Înlocuirea în log-probabilitatea datelor complete generează funcția Q
( Q(beta, sigma^2 mid beta^{
Sau echivalent, folosind notația vectorială pentru vectorul rezidual (y – X beta^{
( Q(beta, sigma^2 mid beta^{
Pasul E se prăbușește pentru a introduce valorile actuale ajustate. Nu este nevoie de integrare.
Etapa de maximizare
Pasul M actualizează parametrii prin maximizarea (Q) în raport cu (beta) (și (sigma^2)).
Se actualizează (beta). Singurul termen din (Q) care depinde de (beta) este suma reziduurilor pătrate. Aplicând regula lanțului la (Q) în raport cu (beta_j) dă
( frac{partial Q}{partial beta_j} = -frac{1}{2sigma^2} sum_{i=1}^n 2!left(y_i – sum_{k=1}^p beta_k x_{ik}right)(-x_{ij}) = frac{^2}=^_{^2}={sigma{1}} x_{ij}!left(y_i – sum_{k=1}^p beta_k x_{ik}right).
Setând acest lucru la zero și înmulțind prin (sigma^2):
( sum_{i=1}^n x_{ij} y_i = sum_{i=1}^n x_{ij} sum_{k=1}^p beta_k x_{ik} = sum_{k=1}^p beta_k underbrace{sum_{i=1}^n x_{ij}^n x_{ij}_{_{ij}_{_{i}}_{_{k}} quad j = 1, ldots, p.
Partea stângă este (j)-a intrare a lui (X^T y), deoarece ((X^T y)_j = sum_i x_{ij} y_i). Partea dreaptă este (j)-a intrare a lui (X^TX beta), deoarece intrarea ((j,k)) a lui (X^TX) este exact (sum_i x_{ij} x_{ik}). Stivuirea tuturor ecuațiilor (p) ((j = 1, ldots, p)) într-o singură ecuație matriceală:
( X^TX beta = X^T y. )
unde (X^TX in mathbb{R}^{p times p}) este o matrice simetrică pozitiv-definită (presupunând că coloanele lui (X) sunt liniar independente) și (X^T y in mathbb{R}^p) este un vector de produse interne între fiecare caracteristică și răspuns. Pentru ilustrare cu (p = 3), acestea arată astfel:
( underbrace{begin{pmatrix} sum x_{i1}^2 și sum x_{i1}x_{i2} și sum x_{i1}x_{i3} \ sum x_{i1}x_{i2} și sum x_{i2}^2 și sum x_{i2}x_{i2} x_{i2} x_{i2} x_{3} sum x_{i2}x_{i3} și sum x_{i3}^2 end{pmatrix}}_{X^TX} begin{pmatrix} beta_1 \ beta_2 \ beta_3 end{pmatrix} = underbrace{begin{pmatrix} sum y_{ii1} sum x_{ii} sum x_{i3} y_i end{pmatrix}}_{X^T y}.
Inversarea (X^TX) dă formula MOL
( beta^{(t+1)} = (X^TX)^{-1} X^T y. )
Rețineți că această soluție este independentă de iterația curentă (beta^{
Se actualizează (sigma^2). Setarea (partial Q / partial sigma^2 = 0):
( sigma^{2(t+1)} = frac{1}{n} sum_{i=1}^n left(y_i – sum_{j=1}^p beta_j^{(t+1)} x_{ij}right)^2 = frac{1}{n}|y – X beta_{1}^{(t+1)}
care este media reziduală pătratică la ponderile actualizate.
Convergenţă
Deoarece pasul M produce maximul global al funcției Q în formă închisă și acel maxim nu depinde de (beta^{
Modelele de regresie Poisson numără datele. În afară de aceasta, modelele de comerț internațional se bazează pe pseudo-probabilitatea maximă Poisson (PPML) și o variabilă continuă, cum ar fi exporturile (sau importurile) nu urmează o distribuție Poisson discretă, motiv pentru care numele PPML și nu PML. Estimatorul PPML este consecvent dacă media condiționată a variabilei de interes este corect specificată. Mai multe despre asta pe pagina The Log of Gravity.
Răspunsul (y_i in {0, 1, 2, ldots}) se presupune că urmează o distribuție Poisson a cărei medie depinde de covariate printr-un log-link:
( y_i sim text{Poisson}(mu_i), quad mu_i = exp!left(sum_{j=1}^p beta_j x_{ij}right) = exp(x_i^T beta). )
Spre deosebire de regresia liniară, nu există o soluție în formă închisă pentru (beta), așa că avem nevoie de o metodă iterativă. Algoritmul EM oferă unul prin introducerea de variabile latente care fac problema de date complete tratabilă.
Construcție latentă-variabilă
Scrieți media Poisson ca (mu_i = exp(x_i^T beta)) și introduceți (m) indicatori binari latenți (z_{i1}, ldots, z_{im}). Acesta este unul pentru fiecare dintre (m) subprocese ipotetice care împreună generează (y_i). Mai exact, împărțiți (mu_i) în (m) părți egale (lambda = mu_i / m) și lăsați
( z_{il} sim text{Bernoulli}(lambda / (1 + lambda)), quad l = 1, ldots, m, )
astfel încât (y_i = sum_{l=1}^m z_{il}) în limita (m to infty).
În practică, formularea standard EM pentru regresia Poisson evită această construcție explicită și tratează în schimb datele complete ca perechea ((y_i, eta_i)), unde (eta_i = x_i^T beta) este predictorul liniar și exploatează proprietățile familiei exponențiale care se aplică log-probabilității Poisson.
Pentru mai multe despre familia exponențială și suficiența statistică, puteți verifica Casella și Berger. Este una dintre cărțile mele preferate (spre deosebire de altele care pun eleganța peste claritate).
Configurare: probabilitatea de înregistrare a datelor complete
Log-probabilitatea Poisson pentru o singură observație este
( log p(y_i mid beta) = y_i log mu_i – mu_i – log(y_i!) = y_i (x_i^T beta) – exp(x_i^T beta) – log(y_i!). )
Însumarea tuturor observațiilor (n) dă probabilitatea logaritării datelor complete (scăderea constantei (sum_i log(y_i!))):
( ell(beta) propto sum_{i=1}^n left( y_i (x_i^T beta) – exp(x_i^T beta) right) = y^TX beta – mathbf{1}^T exp(Xbeta), )
unde (exp(Xbeta)) denotă exponentiația în funcție de element și (mathbf{1}) este un vector de unități.
Pas de așteptare
Spre deosebire de regresia liniară, log-probabilitatea Poisson nu este pătratică în (beta), astfel încât pasul E nu se prăbușește trivial. Abordarea standard este de a construi un surogat pătratic de lucru (funcția Q) la iterația curentă (beta^{
( exp(x_i^T beta) approx exp(eta_i^{
În afară de aceasta, expansiunile Taylor oferă fundația pentru metoda lui Newton. Ambele sunt utilizate intens în industrie, iar un exemplu celebru este Quake’s III Fast Inverse Square Root.
Înlocuirea în (ell(beta)) și păstrarea numai a termenilor care depind de (beta) dă funcția Q
( Q(beta mid beta^{
unde (mu_i^{
Definirea răspunsului de lucru
( tilde{y}_i^{
iar greutatea (mu_i^{
( Q(beta mid beta^{
Acesta este exact un obiectiv ponderat al celor mai mici pătrate. Este aceeași structură ca și log-probabilitatea regresiei liniare cu date complete, dar cu ponderi specifice observației (mu_i^{
Etapa de maximizare
Maximizarea funcției Q cu cele mai mici pătrate ponderate în raport cu (beta) este același calcul ca în pasul M de regresie liniară, dar cu o matrice a greutății diagonale (W^{
Se actualizează (beta). Ecuațiile normale ponderate sunt obținute prin același argument de gradient de intrare ca înainte. Pentru fiecare (j = 1, ldots, p):
( frac{partial Q}{partial beta_j} = sum_{i=1}^n mu_i^{
care sub formă de matrice este
( X^TW^{
Aici (X^TW^{
( beta^{(t+1)} = left(X^TW^{
Aceasta este actualizarea Iteratively Reweighted Least Squares (IRLS), care este algoritmul standard pentru potrivirea modelelor Poisson (și alte GLM). Fiecare pas M este o problemă MCO ponderată cu aceeași structură ca rezultatul regresiei liniare, dar ponderile (W^{
Convergenţă
Deoarece log-probabilitatea Poisson (ell(beta)) este strict concavă în (beta) (hesianul (-X^TWX) este negativ-definit dacă (X) are rang de coloană complet), secvența ({beta^{
Mai multe despre Poisson și alte modele
Verificați McCullagh și Nelder. Acoperă modelele Poisson, Logit și liniare generalizate cu un tratament foarte detaliat. Aceasta este o altă dintre cărțile mele preferate.
