牛顿拉夫森逻辑回归

-2

我做牛顿拉夫森逻辑回归的代码。不幸的是，我尝试了很多数据，没有收敛。有一个错误，我不知道它在哪里。任何人都可以帮助弄清楚是什么问题。牛顿拉夫森逻辑回归

首先数据如下; 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

发现它！你正在更新p基于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钳工不工作（它使用迭代重加权最小二乘法，不NR）：（如果你想查看结果，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)

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

牛顿拉夫森逻辑回归

回答

相关问题