Logistic

為什麼使用牛頓法進行邏輯回歸優化稱為迭代重加權最小二乘法?

  • May 3, 2018

為什麼使用牛頓法進行邏輯回歸優化稱為迭代重加權最小二乘法?

我似乎不清楚,因為邏輯損失和最小二乘損失是完全不同的東西。

總結:GLM 是通過Fisher 評分擬合的,正如 Dimitriy V. Masterov 所指出的那樣,它是 Newton-Raphson 與預期的 Hessian(即,我們使用 Fisher 信息的估計而不是觀察到的信息)。如果我們使用規範鏈接函數,結果是觀察到的 Hessian 等於預期的 Hessian,因此在這種情況下 NR 和 Fisher 評分是相同的。無論哪種方式,我們都會看到,Fisher 評分實際上是擬合加權最小二乘線性模型,並且由此得出的係數估計會收斂* 在邏輯回歸似然的最大值上。除了減少對已解決問題的邏輯回歸擬合之外,我們還可以在最終 WLS 擬合上使用線性回歸診斷來了解我們的邏輯回歸。

我將把重點放在邏輯回歸上,但為了更全面地了解 GLM 中的最大似然性,我推薦本章的第 15.3 節,它通過這一節並在更通用的環境中推導出 IRLS(我認為它來自 John Fox 的Applied回歸分析和廣義線性模型)。

見最後評論


似然函數和評分函數

我們將通過迭代以下形式來擬合我們的 GLM

在哪裡是對數似然和將是對數似然的觀察到的或預期的 Hessian。 我們的鏈接函數是一個函數映射條件均值到我們的線性預測器,所以我們的均值模型是. 讓是將線性預測變量映射到平均值的反向鏈接函數。

對於邏輯回歸,我們有一個獨立觀察的伯努利似然,所以

採取衍生品,


使用規範鏈接

現在讓我們假設我們正在使用規範鏈接功能. 然後所以這意味著這簡化為

所以

此外,還在使用,

然後我們有

並註意這沒有任何在裡面,所以(我們將其視為所以唯一隨機的事情是本身)。因此,當我們在邏輯回歸中使用規範鏈接時,我們已經證明 Fisher 評分等效於 Newton-Raphson。也憑藉 將始終是嚴格負定的,儘管在數值上如果太接近了或者那麼我們可能有權重這可以使負半定,因此在計算上是奇異的。 現在創建工作響應 並註意

總之,這意味著我們可以通過迭代來優化對數似然度

和正是對於加權最小二乘回歸在. 檢查這個R

set.seed(123)
p <- 5
n <- 500
x <- matrix(rnorm(n * p), n, p)
betas <- runif(p, -2, 2)
hc <- function(x) 1 /(1 + exp(-x)) # inverse canonical link
p.true <- hc(x %*% betas)
y <- rbinom(n, 1, p.true)

# fitting with our procedure
my_IRLS_canonical <- function(x, y, b.init, hc, tol=1e-8) {
 change <- Inf
 b.old <- b.init
 while(change > tol) {
   eta <- x %*% b.old  # linear predictor
   y.hat <- hc(eta)
   h.prime_eta <- y.hat * (1 - y.hat)
   z <- (y - y.hat) / h.prime_eta

   b.new <- b.old + lm(z ~ x - 1, weights = h.prime_eta)$coef  # WLS regression
   change <- sqrt(sum((b.new - b.old)^2))
   b.old <- b.new
 }
 b.new
}

my_IRLS_canonical(x, y, rep(1,p), hc)
# x1 x2 x3 x4 x5 
# -1.1149687 2.1897992 1.0271298 0.8702975 -1.2074851

glm(y ~ x - 1, family=binomial())$coef
# x1 x2 x3 x4 x5 
# -1.1149687 2.1897992 1.0271298 0.8702975 -1.2074851 

他們同意。


非規範鏈接函數

現在,如果我們不使用規範鏈接,我們將無法簡化在所以變得更加複雜,因此我們通過使用在我們的費舍爾得分中。

這是怎麼回事:我們已經制定了一般所以 Hessian 將是主要的困難。我們需要

通過期望的線性,我們需要做的就是得到是替換每次出現在我們的模型下,它的平均值是. 因此,總和中的每個術語都將包含一個形式的因子

但要真正進行優化,我們需要估計每個,並且在步 是我們最好的猜測。這意味著這將減少到

這意味著我們將使用和

現在讓

並註意如何在規範鏈接下減少到從上一節。這讓我們寫

除了這是現在而不一定是本身,因此這可能與 Newton-Raphson 不同。對所有人 所以除了數字問題將是負定的。 我們有

所以讓我們的新工作響應成為和, 我們有. 我們一起迭代

所以這仍然是一系列 WLS 回歸,但現在不一定是 Newton-Raphson。


我用這種方式寫出來是為了強調與 Newton-Raphson 的聯繫,但人們經常會考慮更新,以便每個新點本身就是 WLS 解決方案,而不是添加到當前點的 WLS 解決方案. 如果我們想這樣做,我們可以執行以下操作:

因此,如果我們這樣做,您會看到工作響應採用以下形式,但這是同一回事。


讓我們通過使用它對與以前相同的模擬數據執行概率回歸來確認它是否有效(這不是規範鏈接,因此我們需要這種更通用的 IRLS 形式)。

my_IRLS_general <- function(x, y, b.init, h, h.prime, tol=1e-8) {
 change <- Inf
 b.old <- b.init
 while(change > tol) {
   eta <- x %*% b.old  # linear predictor
   y.hat <- h(eta)
   h.prime_eta <- h.prime(eta)
   w_star <- h.prime_eta^2 / (y.hat * (1 - y.hat))
   z_star <- (y - y.hat) / h.prime_eta

   b.new <- b.old + lm(z_star ~ x - 1, weights = w_star)$coef  # WLS

   change <- sqrt(sum((b.new - b.old)^2))
   b.old <- b.new
 }
 b.new
}

# probit inverse link and derivative
h_probit <- function(x) pnorm(x, 0, 1)
h.prime_probit <- function(x) dnorm(x, 0, 1)

my_IRLS_general(x, y, rep(0,p), h_probit, h.prime_probit)
# x1 x2 x3 x4 x5 
# -0.6456508 1.2520266 0.5820856 0.4982678 -0.6768585 

glm(y~x-1, family=binomial(link="probit"))$coef
# x1 x2 x3 x4 x5 
# -0.6456490 1.2520241 0.5820835 0.4982663 -0.6768581 

兩人再次同意。


對收斂的評論

最後,關於收斂性的一些快速評論(我會保持簡短,因為這真的很長而且我不是優化專家)。儘管理論上每個是負定的,糟糕的初始條件仍然可以阻止該算法收斂。在上面的概率示例中,將初始條件更改為會b.init=rep(1,p)導致這種情況,這甚至看起來不像是可疑的初始條件。如果您使用該初始化和這些模擬數據逐步完成 IRLS 過程,那麼在第二次通過循環時,會有一些那一輪正好所以權重變得不確定。如果我們在我給出的算法中使用規範鏈接,我們將永遠不會除以獲得未定義的權重,但如果我們遇到一些情況正在接近或者,例如在完美分離的情況下,那麼我們仍然會得到不收斂,因為梯度消失而我們沒有達到任何目標。

引用自:https://stats.stackexchange.com/questions/344309

comments powered by Disqus