使用广义矩量法(GMM)计算逻辑回归参数

机器算法验证 物流 广义矩
2022-01-26 14:37:50

我想计算与逻辑回归非常相似的回归的系数(实际上是带有另一个系数的逻辑回归: 什么时候可以给出我想过使用 GMM 来计算系数,但我不确定我应该使用什么时刻条件。

A1+e(b0+b1x1+b2x2+),
A

任何人都可以帮助我吗?

谢谢!

2个回答

假设,该模型具有伯努利响应变量A1Yi

Pr(Yi=1)=A1+eXib,

其中(可能还有,取决于它被视为常数还是参数)是拟合系数,是观察数据我假设截距项是通过向数据矩阵添加一个常数值为 1 的变量来处理的。bAXii

瞬间条件是:

E[(YiA1+eXib)Xi]=0.

个观察值,我们将其替换为条件的样本对应项:N

m=1Ni=1N[(YiA1+eXib)Xi]=0

这实际上是通过在所有可能的系数值来解决的(下面我们将使用Nelder-Mead 单纯形法来执行此优化)。mmb

借用关于该主题的优秀 R-bloggers 教程,使用 gmm 包在 R 中实现这一点非常简单。例如,让我们使用 iris 数据集,根据其萼片长度和宽度以及花瓣长度和宽度来预测鸢尾花是否是杂色。在这种情况下,我假设是常数并且等于 1:A

dat <- as.matrix(cbind(data.frame(IsVersicolor = as.numeric(iris$Species == "versicolor"), Intercept=1), iris[,1:4]))
head(dat)
#      IsVersicolor Intercept Sepal.Length Sepal.Width Petal.Length Petal.Width
# [1,]            0         1          5.1         3.5          1.4         0.2
# [2,]            0         1          4.9         3.0          1.4         0.2
# [3,]            0         1          4.7         3.2          1.3         0.2
# [4,]            0         1          4.6         3.1          1.5         0.2
# [5,]            0         1          5.0         3.6          1.4         0.2
# [6,]            0         1          5.4         3.9          1.7         0.4

以下是使用逻辑回归拟合的系数:

summary(glm(IsVersicolor~., data=as.data.frame(dat[,-2]), family="binomial"))
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)    7.3785     2.4993   2.952 0.003155 ** 
# Sepal.Length  -0.2454     0.6496  -0.378 0.705634    
# Sepal.Width   -2.7966     0.7835  -3.569 0.000358 ***
# Petal.Length   1.3136     0.6838   1.921 0.054713 .  
# Petal.Width   -2.7783     1.1731  -2.368 0.017868 *  

我们需要使用 gmm 的主要部分是一个返回矩条件的函数,即每个观察(YiA1+eXib)Xii

moments <- function(b, X) {
  A <- 1
  as.vector(X[,1] - A / (1 + exp(-(X[,-1] %*% cbind(b))))) * X[,-1]
}

我们现在可以数值拟合系数,使用线性回归系数作为一个方便的初始点(如上面链接的教程中所建议的):b

init.coef <- lm(IsVersicolor~., data=as.data.frame(dat[,-2]))$coefficients
library(gmm)
fitted <- gmm(moments, x = dat, t0 = init.coef, type = "iterative", crit = 1e-19,
              wmatrix = "optimal", method = "Nelder-Mead",
              control = list(reltol = 1e-19, maxit = 20000))
fitted
#  (Intercept)  Sepal.Length   Sepal.Width  Petal.Length   Petal.Width  
#      7.37849      -0.24536      -2.79657       1.31364      -2.77834  
# 
# Convergence code =  0 

收敛代码 0 表示过程收敛,参数与逻辑回归返回的参数相同。

快速查看 gmm 包源(提供的函数momentEstim.baseGmm.iterativegmm:::.obj1参数)表明 gmm 包正在最小化mm如上所示。以下等效代码optim直接调用 R 函数,执行我们上面通过调用实现的相同优化gmm

gmm.objective <- function(theta, x, momentFun) {
  avg.moment <- colMeans(momentFun(theta, x))
  sum(avg.moment^2)
}
optim(init.coef, gmm.objective, x=dat, momentFun=moments,
      control = list(reltol = 1e-19, maxit = 20000))$par
#  (Intercept) Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
#    7.3784866   -0.2453567   -2.7965681    1.3136433   -2.7783439 

目前的一种直觉条件是:

E[(Yi11+eXiβ)Xi]

您是否希望您的预测误差为零且与自变量不相关。

一旦你想到ϵi=Yi11+eXiβ作为一个错误术语,规范矩条件类似于您在 OLS 中使用的矩。