

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

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

[code]

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

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


CC BY-SA 4.0
2023-08-07T14:52:28.033
95880
2023-08-16T14:42:28.330
2023-08-08T00:09:42.283
3060
47606
question
r state-machines
ロジスティック回帰の最尤推定におけるヘッセ行列の計算ができない


[R Documentation - Matrix
Diagonals](https://search.r-project.org/R/refmans/base/html/diag.html) の
Details の 4.項によると，対角要素をそれぞれ指定して対角行列を作るには引数をベクトルにする必要があるようです。

[code]

        W <- diag(as.vector(p * (1 - p)))
    
[/code]

御参考ですが，[ロジステック回帰 ->
合格の判定（2次元）](https://oku.edu.mie-u.ac.jp/%7Eokumura/stat/logistic.html)
のデータを使わせていただき，ニュートン法で解（`beta`）を求めてみたところ（記述例は下記），リンク記載の `glm()`
を使う場合と同様の結果が得られました。

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

[code]

    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
    
[/code]


CC BY-SA 4.0
2023-08-16T14:42:28.330
95954
2023-08-16T14:42:28.330
54588
95880
answer