0

ロジスティック回帰の最尤推定におけるヘッセ行列の計算についてご教示ください。

以下、R言語での記述ですが、

hessian <- function(beta, X) {
    p <- exp(X %*% beta) / (1 + exp(X %*% beta))
    W <- diag(p * (1 - p))
    return(t(X) %*% W %*% X) # ヘッセ行列
}

ここで、W <- diag(p * (1 - p)) がスカラーのため、ヘッセ行列が計算できないのですが、どこを修正すれば良いのか分かっておりません。何卒ご教示頂ければ助かります。

1 件の回答 1

0

R Documentation - Matrix Diagonals の Details の 4.項によると,対角要素をそれぞれ指定して対角行列を作るには引数をベクトルにする必要があるようです。

    W <- diag(as.vector(p * (1 - p)))

御参考ですが,ロジステック回帰 -> 合格の判定(2次元) のデータを使わせていただき,ニュートン法で解(beta)を求めてみたところ(記述例は下記),リンク記載の glm() を使う場合と同様の結果が得られました。

なお,pbeta の更新式にも必要なため,求める式(p <- 1 / (1 + exp(-(X %*% beta))))を関数の外に記述し結果を関数に引数で与える形に変更しています。

hessian <- function(p, X) {
    W <- diag(as.vector(p * (1 - p)))
    return(t(X) %*% W %*% X)
}

y  <- c(0, 1, 0, 0, 1, 0, 1, 0, 1, 0)
x1 <- c(3.6, 4.1, 3.7, 4.9, 4.4, 3.6, 4.9, 4.8, 4.1, 4.3)
x2 <- c(60.1, 52.0, 62.5, 60.6, 54.1, 63.6, 68.0, 38.2, 59.5, 47.3)

X <- cbind(rep(1, length(x1)), x1, x2)
beta <- cbind(rep(0, ncol(X)))

for (i in 1:5) {
    p <- 1 / (1 + exp(-(X %*% beta)))
    H <- hessian(p, X)
    beta <- beta - solve(H) %*% t(X) %*% (p - y)
    print(t(beta))
}
##                      x1         x2
## [1,] -8.777639 1.183927 0.05933535
##                      x1         x2
## [1,] -9.435798 1.270301 0.06415955
##                      x1         x2
## [1,] -9.445889 1.271584 0.06424048
##                      x1        x2
## [1,] -9.445892 1.271584 0.0642405
##                      x1        x2
## [1,] -9.445892 1.271584 0.0642405

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。