Am învățat să obținem algoritmul Newton-Raphson din aproximarea seriei Taylor și îl punem în aplicare pentru regresia logistică în R. Vom arăta cum extinderea Taylor de ordinul doi duce la formula de actualizare Newton-Raphson, apoi să comparăm actualizări ale parametrilor individuali față de utilizarea Matricului informațional Full Fisher pentru o convergență mai rapidă.
Obiective
Dovada aproximării seriei Taylor la algoritmul Newton Raphson
Ce este seria Taylor?
O aproximare a seriei Taylor este o tehnică matematică care reprezintă o funcție lină ca o sumă infinită de termeni calculați din derivatele funcției într -un singur punct. Numită după matematicianul Brook Taylor, această metodă exprimă o funcție f (x) în jurul unui punct „a” ca f(a) + f'(a)(x-a) + f''(a)(x-a)²/2! + f'''(a)(x-a)³/3! + ...unde fiecare termen implică instrumente derivate și puteri de (XA). Frumusețea seriei Taylor constă în capacitatea lor de a aproxima funcții complexe folosind termeni polinomiali simpli – cu cât includeți mai mulți termeni, cu atât aproximarea dvs. mai exactă devine într -o anumită rază de convergență în jurul punctului de expansiune. Acest lucru face ca seria Taylor să fie de neprețuit în calcul, fizică și inginerie pentru rezolvarea ecuațiilor diferențiale, analizarea oscilațiilor și efectuarea de calcule numerice în care soluțiile exacte sunt dificil de obținut.
Ce este algoritmul Newton Raphson?
Algoritmul Newton-Raphson este o metodă numerică iterativă folosită pentru a găsi aproximări succesiv mai bune ale rădăcinilor (sau zero) ale unei funcții cu valoare reală. Se bazează pe ideea că, dacă aveți o funcție f (x) și derivatul său f ‘(x), puteți utiliza linia tangentă într -un punct x₀ pentru a găsi o aproximare mai bună a rădăcinii. Formula pentru următoarea aproximare X₁ este dată de x₁ = x₀ - f(x₀)/f'(x₀). Prin repetarea acestui proces, puteți converge la rădăcina reală a funcției. Metoda este deosebit de eficientă pentru funcțiile care sunt continue și diferențiate și converg rapid atunci când ghicirea inițială este apropiată de adevărata rădăcină.
Cum obținem de la seria Taylor la Newton Raphson – Dovadă
Teorema 1: Teorema Taylor
$$ f(x+h) = f(x) + f'(x)h + frac{1}{2}f”(x)h^2 + frac{1}{3!}f”'(x)h^3 + ldots + frac{1}{k!}f^{(k)}(x)h^k + frac {1} {(k+1)!} f^{(k+1)} (w) h^{k+1} $$ Se poate arăta că, pe măsură ce h merge la 0 termenii de ordin superior din teorema noastră Taylor, mergeți la 0 mult mai repede decât H merge la 0.
Dacă ar fi să folosim aproximarea Taylor de ordinul doi, avem: $$ f (x + h) ≈ f (x) + f ‘(x) h + frac {1} {2} f’ ‘(x) h^2 $$, deci avem o f(x) Și vrem să adăugăm o valoare h la x. Extinderea seriei Taylor ne oferă o aproximare a function f(x + h)
Să numim această aproximare g(h)
$$ g (h) = f (x) + f ‘(x) h + frac {1} {2} f’ ‘(x) h^2 $$ Vrem să găsim valoarea din h asta maximizează g(h). Asta înseamnă că putem lua derivatul lui G (h) în raport cu H și să -l setăm egal cu 0.
$$ begin {adună} frac { parțial} { parțial h} g (h) = 0 + f ‘(x) + f’ ‘(x) h \ 0 = f’ (x) + f ” (x) h \ -f ‘(x) = f’ ‘(x) h \ h = – frac {f noi valoarea h asta maximizează g(h).
Dacă vrem să estimăm noua valoare a xputem adăuga h la x: $$ x_ {new} = x_ {old} + h \ = x_ {old} – frac {f ‘(x_ {old})} {f’ ‘(x_ {old})} $$ et viola! De la seria Taylor la alogoritmul Newton Raphson! ❤️
Scufundați -vă mai adânc
Acum, să continuăm călătoria noastră de regresie logistică prin implementarea algoritmului Newton-Raphson pentru a estima coeficienții folosind MLE.
Anterior, am stabilit funcția de probabilitate pentru regresia logistică și am derivat funcția de probabilitate de log pentru beta0 (interceptare) și beta1 (coeficient pentru 1 predictor).
$$ ln L(boldsymbol{beta_0, beta_1}) = sum_{i=1}^{n} left( y_i ln(p_i) + (1-y_i) ln(1-p_i) right) $$ For beta0, to find the derivate of the log-likelihood function, we need to take the first derivative of the Funcția de lungă durată de jurnal în raport cu Beta0, ca atare.
$$ frac { partial ln l ( Boldsymbol { beta_0, beta_1})} { partial beta_0} = sum_ {i = 1}^{n} left (y_i-p_i dreapta) $ $ pentru beta1, trebuie să luați primul derivativ al funcției log-likeli $$ frac { partial ln l ( Boldsymbol { beta_0, beta_1})} { partial beta_1} = sum_ {i = 1}^{n} left (y_i – p_i dreapta) x_ {1i} $ $ vor fi noi f'(x)
Al doilea derivat al funcției de log -likelihood cu privire la beta0 este dat de: $$ frac { parțial^2 ln l ( Boldsymbol { beta_0, beta_1})} { parțial beta_0^2} = – sum_ {i = 1}^{n} p_i (1 -p_) $ $ $ pentru Al doilea derivat al funcției de probabilitate a jurnalului în raport cu Beta1, așa cum este așa. $$ frac { partial^2 ln l ( Boldsymbol { beta_0, beta_1})} { parțial beta_1^2} = – sum_ {i = 1}^{n} p_i (1 -p_i) x_ {1i}^2 $ $
Acestea vor fi ale noastre f''(x)
Să punem totul împreună
Acum putem folosi algoritmul Newton-Raphson pentru a estima coeficienții modelului de regresie logistică. Algoritmul va actualiza iterativ coeficienții până la obținerea convergenței.
Beta0
$$ begin{gather} beta_{new} = beta_{old} + h \ = beta_{old} – frac{f'(beta_{old})}{f”(beta_{old})} \ = beta_{old} – (frac{sum_{i=1}^{n} left (y_i-p_i dreapta)} {- sum_ {i = 1}^{n} p_i (1-p_i)}) \ = beta_ {old} + frac { sum_ {i = 1}^{n} left (y_i-p_i dreapta)} { sum_ {i = 1}^{n}}} { sum_ {i = 1}^{n}}} { sum_ {i = 1}^{n}}} { sum_ {i = 1}^{n} p_ (1 – p_i)} end {adună} $$
Beta1
$$ begin{gather} beta_{new} = beta_{old} + h \ = beta_{old} – frac{f'(beta_{old})}{f”(beta_{old})} \ = beta_{old} – (frac{sum_{i=1}^{n} left (y_i-p_i dreapta) x_ {1i}} {- sum_ {i = 1}^{n} p_i (1-p_i) x_ {1i}^2}) \ = beta_ {vechi} + frac { sum_ {i = 1}^{n} tet x_ {1i}} { sum_ {i = 1}^{n} p_i (1 – p_i) x_ {1i}^2} end {adun} $$ wow, uimitor! Am făcut -o! Acum, să punem aceste formule în cod și să vedem cum funcționează.
Să codificăm
library(tidyverse)
set.seed(100)
n <- 100
x <- rbinom(n,1,0.5)
y <- rbinom(n,1,plogis(-2+2*x))
beta <- c(0,0)
iter <- 100
history <- matrix(0,nrow = iter, ncol = 2)
tolerance <- 10^-8
for (i in 1:iter) {
z <- beta(1) + beta(2)*x
p <- 1 / (1+exp(-z))
d_b0 <- sum(y-p)
d_b1 <- sum((y-p)*x)
d2_b0 <- -sum(p*(1-p))
d2_b1 <- -sum(p*(1-p)*x^2)
beta_new <- beta - c(d_b0/d2_b0, d_b1/d2_b1)
history(i,) <- beta_new
if (abs(beta_new(1) - beta(1)) < tolerance && abs(beta_new(2) - beta(2)) < tolerance) {
cat(paste0("converged on iter ", i,"nbeta0: ",beta_new(1),"nbeta1: ",beta_new(2)))
break
}
beta <- beta_new
if (i==iter) {
cat(paste0("converged on iter ", i,"nbeta0: ",beta_new(1),"nbeta1: ",beta_new(2)))
}
}
## converged on iter 100
## beta0: -2.75149081264213
## beta1: 2.99265005420773
Bine! Am implementat cu succes algoritmul Newton-Raphson pentru a estima coeficienții modelului de regresie logistică. Algoritmul actualizează iterativ coeficienții până la obținerea convergenței și putem vedea estimările finale ale beta0 și beta1.
Este ciudat că a fost nevoie de atât de multe iterații pentru a converge decât GLM, deoarece setul maxim de iter pentru GLM este de 25 de ani.
glm.control
## function (epsilon = 1e-08, maxit = 25, trace = FALSE)
## {
## if (!is.numeric(epsilon) || epsilon <= 0)
## stop("value of 'epsilon' must be > 0")
## if (!is.numeric(maxit) || maxit <= 0)
## stop("maximum number of iterations must be > 0")
## list(epsilon = epsilon, maxit = maxit, trace = trace)
## }
##
##
Comparați cu GLM
(summary(glm(y~x,family = binomial(link = "logit")))) ## ## Call: ## glm(formula = y ~ x, family = binomial(link = "logit")) ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.7515 0.5955 -4.621 3.83e-06 *** ## x 2.9927 0.6601 4.533 5.80e-06 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 123.82 on 99 degrees of freedom ## Residual deviance: 91.29 on 98 degrees of freedom ## AIC: 95.29 ## ## Number of Fisher Scoring iterations: 5
OK, estimările punctului arată foarte aproape unul de celălalt. Dar de ce este iterația noastră cu atât mai mult decât GLM? Acest lucru se datorează faptului că nu folosim matricea completă a informațiilor Fisher pentru a ne actualiza coeficienții.
Matricea informațională Fisher
După cum s-a menționat anterior, matricea informațională Fisher este o matrice pătrată care conține derivate parțiale de ordinul doi ale funcției de probabilitate de jurnal în raport cu parametrii. Oferă informații despre curbura funcției de probabilitate a jurnalului și este utilizat pentru a estima matricea de varianță-covarianță a estimărilor de probabilitate maximă. Nu numai că, folosind în ansamblu, optimizează în continuare convergența algoritmului.
$$ begin {adună} i ( Boldsymbol { beta}) = begin {bmatrix} frac { parțial^2 ln l ( Boldsymbol { beta})} { parțial beta_0^2} & frac { parțial^2 ln L ( Boldsymbol { beta})} { parțial beta_0 parțial beta_1} \ frac { parțial^2 ln l ( Boldsymbol { beta})} { parțial beta_1 parțial beta_0} & frac { parțial^2 ln L ( Boldsymbol { beta})} { partial beta_1^2} end {bmatrix} \ = begin {bmatrix} – sum_ {i = 1}^{n} p_i (1 – p_i) & – sum_ {i = 1}^{n} p_i (1 – p_) x_ {1i} \ – sum_ {i = 1}^{n} p_i (1 – p_i) x_ {1i} & – sum_ {i = 1}^{n} p_i (1 – p_i) x_ {1i}^2 end {bmatrix} end {adunat} $ $ $ $
Să punem totul împreună
Acum putem folosi matricea de informare Fisher pentru a actualiza coeficienții modelului de regresie logistică. Algoritmul va actualiza iterativ coeficienții până la obținerea convergenței.
Să extindem formulele noastre anterioare pentru a include matricea de informare Fisher și vectorul scorului (aka gradient al vectorilor) ca așa.
$$ beta_ {new} = beta_ {old} + h \ begin {bmatrix} beta_ {0_ text {new}} \ beta_ {1_ text {new}} end {bmatrix} = begin {bmatrix} beta_ {0_ text {vechi {bmatrix} beta_ {0 beta_ {1_ text {old}} end {bmatrix} – i ( Boldsymbol { beta_ {Old}})^{ – 1} cdot nabla f ( beta_ {Old}) \ = begin {BMatrix} beta_ {0_ text {Old}} beta_ {1_ text {old}} end {bmatrix} – begin {bmatrix} frac { parțial^2 ln l ( Boldsymbol { beta})} { parțial beta_0^2} & frac { parțial^2 ln L ( Boldsymbol { beta})} { parțial beta_0 parțial beta_1} \ frac { parțial^2 ln l ( Boldsymbol { beta})} { parțial beta_1 parțial beta_0} & frac { parțial^2 ln L ( Boldsymbol { beta})} { parțial beta_1^2} end {bmatrix}^{-1} cdot begin {bmatrix} frac { parțial ln l ( Boldsymbol { beta}) {{ parțial beta_0} beta}) { parțial ln L ( Boldsymbol { beta})} { parțial beta_1} end {bmatrix} $$ Să încercăm să facem un pic mai puțin messy $$ text {let a} = frac { parțial^2 ln l ( Boldsymbol { beta})} { parțial beta_0^2 ; text{d} = frac{partial^2 ln L(boldsymbol{beta})}{partial beta_1^2} \ begin{bmatrix} beta_{0_text{new}} \ beta_{1_text{new}} end{bmatrix} = begin{bmatrix} beta_{0_text{old}} \ beta_{1_text{old}} end{bmatrix} – begin{bmatrix} a & b \ c & d end{bmatrix}^{-1} cdot begin{bmatrix} frac{partial ln L(boldsymbol{beta})}{partial beta_0} \ frac{partial ln L(boldsymbol{beta})}{partial beta_1} end{bmatrix} \ = begin{bmatrix} beta_{0_text{old}} \ beta_ {1_ text {vechi}} end {bmatrix} – frac {1} {ad – bc} begin {bmatrix} d & -b \ -c & a end {bmatrix} cdot begin {bmatrix} frac { parțial lnot begin {bmatrix} frac { parțial L ( Boldsymbol { beta})} { partial beta_0} \ frac { partial ln l ( Boldsymbol { beta})} { partial beta_1} end {bmatrix} text {BMatrix beta_ {1_ text {Old}} end {bmatrix} – begin {bmatrix} frac {d} {ad – bc} & frac {-b} {ad – bc} \ frac {-c} {ad – bc} & frac {a} {ad – bc}} & frac {a} {ad – Bc} end {bmatrix} cdot begin {bmatrix} frac { parțial ln l ( Boldsymbol { beta})} { parțial beta_0} \ frac { parțial ln l ( beta_1 { beta})} { parțial beta_1 { beta})} { parțial end {bmatrix} $$ După cum putem vedea aici, avem o înmulțire matricială a inversului matricei informaționale Fisher și a vectorului scorului. În exemplul nostru anterior simplu, am inclus doar elementele individuale diagonale ale matricei informaționale Fisher. Prin includerea elementelor off-diagonale, putem obține o estimare mai exactă a coeficienților și probabil o convergență mai rapidă ca în GLM. Să codificăm și să vedem dacă este adevărat!
Iterarea informației Fisher
beta <- c(0,0)
iter <- 25
history <- matrix(0,nrow = iter, ncol = 2)
tolerance <- 10^-8
for (i in 1:iter) {
z <- beta(1) + beta(2)*x
p <- 1 / (1+exp(-z))
d_b0 <- sum(y-p)
d_b1 <- sum((y-p)*x)
score_vec <- c(d_b0, d_b1)
i_11 <- sum(p*(1-p))
i_10 <- sum(x*p*(1-p))
i_01 <- i_10
i_22 <- sum(x^2*p*(1-p))
i_mat <- matrix(c(i_11,i_01,i_10,i_22),nrow = 2, ncol = 2)
i_mat_inv <- solve(i_mat)
beta_new <- beta + i_mat_inv %*% score_vec
history(i,) <- beta_new
## se
se_beta0 <- sqrt(diag(i_mat_inv)(1))
se_beta1 <- sqrt(diag(i_mat_inv)(2))
if (abs(beta_new(1) - beta(1)) < tolerance && abs(beta_new(2) - beta(2)) < tolerance) {
cat(paste0("converged on iter ", i,"nbeta0: ",beta_new(1)," (",se_beta0,") ","nbeta1: ",beta_new(2)," (",se_beta1,") "))
break
}
beta <- beta_new
if (i==iter) {
cat(paste0("converged on iter ", i,"nbeta0: ",beta_new(1)," (",se_beta0,") ","nbeta1: ",beta_new(2)," (",se_beta1,") "))
}
}
## converged on iter 7
## beta0: -2.75153531304195 (0.595491334175413)
## beta1: 2.99269736985884 (0.660135410538508)
Uau, nu prea shabby! Estimările punctuale și SE sunt foarte aproape de GLM. Iterațiile au fost cu siguranță mai scurte decât exemplul anterior.
Acolo îl avem! Este fascinant să privim sub capotă cum funcționează toate acestea. Algoritmul Newton-Raphson este un instrument puternic pentru estimarea coeficienților în regresia logistică și înțelegerea conexiunii sale la aproximarea seriei Taylor oferă o perspectivă valoroasă asupra proprietăților sale de convergență.
Lecții învățate
- A învățat înțelegerea de suprafață a ceea ce face aproximarea seriei Taylor.
- A învățat algoritmul Newton-Raphson. O astfel de formulă elegantă!
- Dovada algoritmului Newton-Raphson folosind aproximarea seriei Taylor
- Cum includerea matricei informaționale Fisher poate îmbunătăți convergența, anterior am crezut întotdeauna că coeficienții sunt independenți unul de celălalt și sunt necesare doar elementele diagonale individuale. Dar se dovedește că elementele off-diagonale sunt de asemenea importante pentru convergență.
- a învățat
(nabla)Simbolul se numește simbolul NABLA, care este utilizat pentru a denota gradientul unei funcții. - învăţat
glm.controlși modul în care Maxiterul implicit GLM este 25 și Epsilon (pragul de toleranță) este 10^-8.
Dacă vă place acest articol:
