Использование обобщенного метода моментов (GMM) для расчета параметра логистической регрессии

Я хочу вычислить коэффициенты для регрессии, которая очень похожа на логистическую регрессию (На самом деле логистическая регрессия с другим коэффициентом: когда может быть дано). Я думал об использовании GMM для вычисления коэффициентов, но я не уверен, какие условия момента я должен использовать.

\frac{A}{1 + e^{- (b_{0} + b_{1} x_{1} + b_{2} x_{2} + \dots)}},

$\frac{A}{1 + e^{- (b_0 + b_1 x_1 + b_2 x_2 + \ldots)}},$

A

$A$

Кто-нибудь может мне помочь с этим?

Благодарность!

logistic generalized-moments user5497
источник

Когда вы говорите «

может быть дано», вы имеете в виду, что он указан пользователем или оценивается моделью?

A

$A$

Макро

так или иначе. Я могу поставить его как вход (например, А = 0,25) или быть одним из коэффициентов, которые будут найдены

user5497

Отличается ли он от субъекта к субъекту (т.е. это данные) или это фиксированная константа во всех наблюдениях?

Макро

исправлено во всех наблюдениях (например, b0, b1, ...)

user5497

Почему бы не использовать максимальную вероятность вместо GMM?

Макро

Ответы:

Предполагая , эта модель имеет переменную отклика Бернулли с $A\leq 1$ $Y_i$

P r (Y_{i} = 1) = \frac{A}{1 + e^{- X_{i}^{'} b}},

$Pr(Y_i = 1) = \frac{A}{1+e^{-X_i'b}},$

где (и, возможно, , в зависимости от того, рассматривается ли это как константа или параметр) - это подогнанные коэффициенты, а - данные для наблюдения . Я предполагаю, что термин перехват обрабатывается путем добавления переменной с постоянным значением 1 в матрицу данных. $b$ $A$ $X_i$ $i$

Моментальные условия:

\begin{aligned} E [(Y_{i} - \frac{A}{1 + e^{- X_{i}^{'} b}}) X_{i}] & = 0. \end{aligned}

$\begin{align*} \mathbb{E}\bigg[\bigg(Y_i-\frac{A}{1+e^{-X_i'b}}\bigg)X_i\bigg] &= 0. \end{align*}$

$N$

m = \frac{1}{N} \sum_{i = 1}^{N} [(Y_{i} - \frac{A}{1 + e^{- X_{i}^{'} b}}) X_{i}] = 0

$m = \frac{1}{N}\sum_{i=1}^N \bigg[\bigg(Y_i-\frac{A}{1+e^{-X_i'b}}\bigg)X_i\bigg] = 0$

$m'm$ $b$

$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 *

$(Y_i-\frac{A}{1+e^{-X_i'b}})X_i$ $i$

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 указывает на сходящуюся процедуру, а параметры идентичны тем, которые возвращаются логистической регрессией.

momentEstim.baseGmm.iterativegmm:::.obj1 $m'm$ optimgmm

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

josliber
источник