数值方法求解逻辑回归

Cloud2016

模拟一个数据集，数据是根据逻辑回归模型生成的

set.seed(2023)
n <- 2500
k <- 10
X <- matrix(rnorm(n * k), ncol = k)
y <- rbinom(n, size = 1, prob = plogis(1 + 3 * X[,1] - 2 * X[,2]))

现在想根据这组数据，求解逻辑回归模型的系数，不要用现成函数 glm.fit()，把它当作一个无约束优化问题来求解。

log_logit_lik <- function(beta) {
  -sum(y * log(plogis(cbind(1, X) %*% beta)) - (1 - y) * log(1 - plogis(cbind(1, X) %*% beta)))
}

当我用 optim() 来求解发现，结果距离真值差距非常大，感觉负对数似然函数写错了，但是又不知道错哪了？求助

# Nelder-Mead
optim(
  par = rep(1, 11),   # 初始值
  fn = log_logit_lik, # 目标函数
  method = "Nelder-Mead"
)

$par
 [1] 11.43868201 -3.91046434  4.99817589 -1.29366216 -1.25622057 -0.81202657
 [7]  1.83220909  0.09821513  5.02672298 -1.31117426  1.60862047

$value
[1] -16007.3

$counts
function gradient 
     502       NA 

$convergence
[1] 1

$message
NULL

逻辑回归模型的参考书 https://www.stat.umn.edu/geyer/5931/mle/glm.pdf

没明白 glm() 里面是咋算的，好像也没什么材料说计算过程，求助一下。

yihui

Cloud2016 似然函数写错了，1 - y 前面的减号应该是加号。要是不对着数学公式直接抄的话，代码可以写得简练、高效（避免重复计算矩阵乘法 X * beta），这样就更容易看出问题来了：

log_logit_lik <- function(beta) {
  p <- plogis(cbind(1, X) %*% beta)
  -sum(y * log(p) + (1 - y) * log(1 - p))
}

而且 optim() 的结果都没有收敛，convergence = 1，没收敛的情况下还无法下结论。1 表示到达了最大迭代次数还不收敛，所以应该增加迭代次数。我试了一下，用模拟退火以及一万次迭代的结果比较能接近真实参数（c(1, 3, -2, 0, 0, 0, ...)）。

optim(
    par = rep(1, 11),   # 初始值
    fn = log_logit_lik, # 目标函数
    method = "SANN", 
    control = list(maxit = 10000)
)

至于 glm() 是怎么算的，当年我还没毕业就已经好心地提前还给老师了。模糊记得有个什么迭代加权最小二乘（IRLS），但现在打死我也想不起来那是啥玩意儿了。实在要搞清楚的话，你可以看 glm() 的源代码。

Cloud2016 逻辑回归模型的参考书 https://www.stat.umn.edu/geyer/5931/mle/glm.pdf

又及：说起这个 Charles J. Geyer 教授，他在统计计算方面是很强，也可谓是可重复性研究的先驱和倡导者之一：http://users.stat.umn.edu/~geyer/Sweave/ 而且他很多年前捣鼓的 Rweb 也是极其有开创性：https://rweb.webapps.cla.umn.edu/Rweb/ 简直可谓 RStudio 和 OpenCPU 的前身。当年看得我眼前一亮，还写了一封信去膜拜和咨询几个技术问题，结果没人鸟我，可能那时候人微言轻、英语又菜吧，哈哈。但直到现在我还是有点膜拜这位老爷爷的。

Cloud2016

yihui 我其实有公式，可能是狗眼盯着看久了，硬是没看到负号问题。

Cloud2016

yihui 似然函数正确后，可以计算梯度啦，然后利用上梯度信息，迭代效率就快了。

# 梯度函数
log_logit_lik_grad <- function(beta) {
  p <- plogis(cbind(1, X) %*% beta)
  - t((y / p - (1 - y) / (1 - p)) * p * (1 - p))  %*% cbind(1, X) 
}

optim(
  par = rep(1, 11),   # 初始值
  fn = log_logit_lik, # 目标函数
  gr = log_logit_lik_grad, # 目标函数的梯度
  method = "L-BFGS-B"
)

#> $par
#>  [1]  1.00802641  3.11296713 -2.00955313  0.05855394 -0.02650585  0.01330428
#>  [7]  0.02171815  0.10213455 -0.02949774 -0.08633384  0.08098888
#> 
#> $value
#> [1] 750.9724
#> 
#> $counts
#> function gradient 
#>       13       13 
#> 
#> $convergence
#> [1] 0
#> 
#> $message
#> [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

从负对数似然函数的值来看，更准了。

yihui

Cloud2016 盯着看久了，硬是没看到负号问题

早睡早起。这都几点了。

Cloud2016 狗眼

哈哈，很好。

Cloud2016 利用上梯度信息，迭代效率就快了

人脑一小步，电脑一大步。很好。

nan.xiao

如果楼主也像我一样是懒人想直接抄公式的话…… 把问题在计算图上定义好，直接 SGD 求解就可以了，有自动求导，不需要再去手动推导梯度。Cross entropy loss 单独实现这样比较清楚而且可以复用。把每次迭代的误差项画个图也可以看出实现或者收敛是否有问题。

set.seed(2023)
n <- 2500
k <- 10
X <- matrix(rnorm(n * k), ncol = k)
y <- rbinom(n, size = 1, prob = plogis(1 + 3 * X[, 1] - 2 * X[, 2]))

remotes::install_version("cgraph", "6.0.1")
remotes::install_github("nanxstats/logreg")

fit <- logreg::fit_logistic(X, y, n_epochs = 5000)

logreg::plot_error(fit)
logreg::plot_coef(fit)

cgraph::cg_graph_get(fit$graph, "beta")$value
cgraph::cg_graph_get(fit$graph, "alpha")$value

Cloud2016

nan.xiao 这个工具太厉害了，杀鸡就不用牛刀了。如果能手（或者借助符号计算软件）推出来梯度，输入梯度信息，应该还是会比自动微分的方法快吧？

nan.xiao

Cloud2016 如果梯度简单到能手推那当然好了。自动求导框架解决的主要是不方便手推时如何节约人力的自动化问题吧。至于性能就不好评价了，感觉不同框架的实现方差有点大，大概也和具体的优化算法怎么利用梯度信息有关。