你能簡單直觀地解釋一下 IRLS 方法來找到 GLM 的 MLE 嗎？

背景：

我正在嘗試關注普林斯頓對 GLM 的 MLE 估計的評論。

我了解 MLE 估計的基礎知識：likelihood、score、 觀察和預期Fisher information以及Fisher scoring技術。而且我知道如何用 MLE 估計證明簡單的線性回歸。

---

問題：

我什至無法理解這種方法的第一行:(

背後的直覺是什麼工作變量定義為：

為什麼使用它們而不是估計?

它們之間的關係是什麼response/link function？和

如果有人有一個簡單的解釋，或者可以將我引導到一個更基本的文本，我將不勝感激。

幾年前，我為我的學生寫了一篇關於這個的論文（西班牙語），所以我可以嘗試在這裡重寫這些解釋。我將通過一系列增加複雜性的示例來研究 IRLS（迭代重加權最小二乘法）。對於第一個示例，我們需要位置規模族的概念。讓 f0 在某種意義上是一個以零為中心的密度函數。我們可以通過定義來構造一個密度族 f(x)=f(x;μ,σ)=1σf0(x−μσ)

在哪裡 σ>0 是一個尺度參數並且 μ 是位置參數。在測量誤差模型中，通常將誤差項建模為正態分佈，我們可以使用上面構建的位置尺度族來代替該正態分佈。什麼時候 f0 是標準正態分佈，上面的構造給出了 N(μ,σ) 家庭。

現在我們將在一些簡單的例子中使用 IRLS。首先，我們將在模型中找到 ML（最大似然）估計量 Y1,Y2,…,Yni.i.d

隨著密度 f(y)=1π11+(y−μ)2,y∈R,

柯西分佈位置族 μ （所以這是一個位置族）。但首先是一些符號。的加權最小二乘估計量 μ 是（誰）給的 μ∗=∑ni=1wiyi∑ni=1wi.

在哪裡 wi 是一些權重。我們將看到 ML 估計器 μ 可以用相同的形式表示，用 wi 殘差的一些函數 ϵi=yi−ˆμ.

似然函數由下式給出 L(y;μ)=(1π)nn∏i=111+(yi−μ)2

並且對數似然函數由下式給出 l(y)=−nlog(π)−n∑i=1log(1+(yi−μ)2).

它的導數關於 μ 是 ∂l(y)∂μ=0−∑∂∂μlog(1+(yi−μ)2) =−∑2(yi−μ)1+(yi−μ)2⋅(−1) =∑2ϵi1+ϵ2i

在哪裡 ϵi=yi−μ . 寫 f0(ϵ)=1π11+ϵ2 和 f′0(ϵ)=1π−1⋅2ϵ(1+ϵ2)2 ，我們得到 f′0(ϵ)f0(ϵ)=−1⋅2ϵ(1+ϵ2)211+ϵ2=−2ϵ1+ϵ2.

我們發現 ∂l(y)∂μ=−∑f′0(ϵi)f0(ϵi) =−∑f′0(ϵi)f0(ϵi)⋅(−1ϵi)⋅(−ϵi) =∑wiϵi

我們使用定義的地方 wi=f′0(ϵi)f0(ϵi)⋅(−1ϵi)=−2ϵi1+ϵ2i⋅(−1ϵi)=21+ϵ2i.

記住那個 ϵi=yi−μ 我們得到方程 ∑wiyi=μ∑wi,

這是IRLS的估計方程。注意

1. 權重 wi 總是積極的。
2. 如果殘差很大，我們對相應的觀察給予較少的權重。

為了在實踐中計算 ML 估計量，我們需要一個起始值 ˆμ(0) ，例如，我們可以使用中位數。使用這個值我們計算殘差 ϵ(0)i=yi−ˆμ(0)

和權重 w(0)i=21+ϵ(0)i.

的新價值 ˆμ 是（誰）給的 ˆμ(1)=∑w(0)iyi∑w(0)i.

以這種方式繼續我們定義 ϵ(j)i=yi−ˆμ(j)

和 w(j)i=21+ϵ(j)i.

通過時的估計值 j+1 算法變為 ˆμ(j+1)=∑w(j)iyi∑w(j)i.

繼續直到序列 ˆμ(0),ˆμ(1),…,ˆμ(j),…

收斂。

現在我們用更一般的位置和規模族來研究這個過程， f(y)=1σf0(y−μσ) , 細節較少。讓 Y1,Y2,…,Yn 與上面的密度無關。也定義 ϵi=yi−μσ . 對數似然函數是 l(y)=−n2log(σ2)+∑log(f0(yi−μσ)).

寫作 ν=σ2 ， 注意 ∂ϵi∂μ=−1σ

和 ∂ϵi∂ν=(yi−μ)(1√ν)′=(yi−μ)⋅−12σ3.

計算對數似然導數 ∂l(y)∂μ=∑f′0(ϵi)f0(ϵi)⋅∂ϵi∂μ=∑f′0(ϵi)f0(ϵi)⋅(−1σ)=−1σ∑f′o(ϵi)f0(ϵi)⋅(−1ϵi)(−ϵi)=1σ∑wiϵi

並將其等於零給出與第一個示例相同的估計方程。然後尋找估算器 σ2 ： ∂l(y)∂ν=−n21ν+∑f′0(ϵi)f0(ϵi)⋅∂ϵi∂ν =−n21ν+∑f′0(ϵi)f0(ϵi)⋅(−(yi−μ)2σ3) =−n21ν−121σ2∑f′0(ϵi)f0(ϵi)⋅ϵi =−n21ν−121ν∑f′0(ϵi)f0(ϵi)⋅(−1ϵi)(−ϵi)⋅ϵi =−n21ν+121ν∑wiϵ2i!=0.

導致估計器 ^σ2=1n∑wi(yi−ˆμ)2.

在這種情況下也可以使用上面的迭代算法。

在下文中，我們給出了一個使用 R 的數值示例，用於雙指數模型（具有已知比例）和數據y <- c(-5,-1,0,1,5)。對於該數據，ML 估計器的真實值為 0。初始值為mu <- 0.5。該算法的一次通過是

     iterest <- function(y, mu) {
                  w <- 1/abs(y - mu)
                  weighted.mean(y, w)
                  }

使用此功能，您可以嘗試“手動”進行迭代然後迭代算法可以通過

   mu_0 <- 0.5
   repeat {mu <- iterest(y, mu_0)
           if (abs(mu_0 - mu) < 0.000001) break
           mu_0 <- mu }

練習：如果模型是 tk 具有尺度參數的分佈 σ 顯示迭代由權重給出 wi=k+1k+ϵ2i.

練習：如果密度是邏輯的，則顯示權重由下式給出 w(ϵ)=1−eϵ1+eϵ⋅−1ϵ.

暫時我會把它留在這裡，我將繼續這篇文章。

引用自：https://stats.stackexchange.com/questions/236676