로지스틱 회귀에 대한 Newton Raphson

-2

로지스틱 회귀에 대해 Newton Raphson을 코딩했습니다. 불행히도 나는 많은 데이터를 수렴하지 않았다. 내가 어디 있는지 모르겠다. 누구든지 문제를 파악하는 데 도움을 줄 수 있습니다.로지스틱 회귀에 대한 Newton Raphson

먼저 데이터는 다음과 같습니다. y는 응답 (0,1)을 나타내며, Z는 탐색 변수 인 115 * 30 행렬입니다. 30 가지 매개 변수를 추정해야합니다.

y = c(rep(0,60),rep(1,55)) 
X = sample(c(0,1),size=3450,replace=T) 
Z = t(matrix(X,ncol=115)) 
#The code is ; 
B  = matrix(rep(0,30*10),ncol=10) 
B[,1] = matrix(rep(0,30),ncol=1) 
for(i in 2 : 10){ 
    print(i) 
    p  <- exp(Z %*%as.matrix(B[,i]))/(1 + exp(Z %*% as.matrix(B[,i]))) 
    v.2  <- diag(as.vector(1 * p*(1-p))) 
    score.2 <- t(Z) %*% (y - p) # score function 
    increm <- solve(t(Z) %*% v.2 %*% Z) 
    B[,i] = as.matrix(B[,i-1])+increm%*%score.2 
    if(B[,i]-B[i-1]==matrix(rep(0.0001,30),ncol=1)){ 
    return(B) 
    } 
}

출처

2014-11-08 Neal

발견! 을 B[,i]으로 업데이트하면 B[,i-1]을 사용해야합니다.

답을 찾던 중 코드를 정리하고 함수에 결과를 포함 시켰습니다. R에 내장 된 glm이 작동하는 것 같습니다 (아래 참조). 한 노트는이 방법이 불안정 할 가능성이 있다는 것입니다 : ... (30 개) 예측 만 115 바이너리 응답과 이진 모델을 피팅하고, 어떤 처벌 또는 수축하지 않고, 매우 낙관적이다

set.seed(101) 
n.obs <- 115 
n.zero <- 60 
n.pred <- 30 
y <- c(rep(0,n.zero),rep(1,n.obs-n.zero)) 
X <- sample(c(0,1),size=n.pred*n.obs,replace=TRUE) 
Z <- t(matrix(X,ncol=n.obs))

R의 내장 glm (. 당신이 library("arm"); coefplot(g1)을 결과를 보려는 경우)

g1 <- glm(y~.-1,data.frame(y,Z),family="binomial")

## B_{m+1} = B_m + (X^T V_m X)^{-1} X^T (Y-P_m)

: 벤치 일 (그것은 반복적으로 최소 제곱하지 NR을 reweighted 사용합니다) 않습니다

NRfit 기능 :

NRfit <- function(y,X,start,n.iter=100,tol=1e-4,verbose=TRUE) { 
    ## used X rather than Z just because it's more standard notation 
    n.pred <- ncol(X) 
    B <- matrix(NA,ncol=n.iter, 
       nrow=n.pred) 
    B[,1] <- start 
    for (i in 2:n.iter) { 
     if (verbose) cat(i,"\n") 
     p <- plogis(X %*% B[,i-1]) 
     v.2 <- diag(c(p*(1-p))) 
     score.2 <- t(X) %*% (y - p) # score function 
     increm <- solve(t(X) %*% v.2 %*% X) 
     B[,i] <- B[,i-1]+increm%*%score.2 
     if (all(abs(B[,i]-B[,i-1]) < tol)) return(B) 
    } 
    B 
} 

matplot(res1 <- t(NRfit(y,Z,start=coef(g1)))) 
matplot(res2 <- t(NRfit(y,Z,start=rep(0,ncol(Z))))) 
all.equal(res2[6,],unname(coef(g1))) ## TRUE

출처

2014-11-08 21:04:54

로지스틱 회귀에 대한 Newton Raphson

답변

관련 문제