Ридж наказал GLM, используя увеличение строки?

Я читал, что регрессия гребня может быть достигнута простым добавлением строк данных в исходную матрицу данных, где каждая строка создается с использованием 0 для зависимых переменных и квадратного корня из $k$ или нуля для независимых переменных. Затем добавляется одна дополнительная строка для каждой независимой переменной.

Мне было интересно, можно ли получить доказательства для всех случаев, в том числе для логистической регрессии или других GLM.

Ридж регрессии минимизирует .$\sum_{i=1}^n (y_i-x_i^T\beta)^2+\lambda\sum_{j=1}^p\beta_j^2$

(Часто требуется постоянная, но не сокращенная. В этом случае она включается в и предикторы - но если вы не хотите сокращать ее, у вас нет соответствующей строки для псевдонаблюдения. Или если вы хотите , чтобы уменьшить его, вы делаете иметь строку для него. Я напишу его , как будто это не учитывается в р , и не сморщенные, так как это более сложный случай. другой случай тривиальное изменение от этого. )$\beta$$p$

Мы можем записать второе слагаемое как псевдонаблюдения, если мы можем записать каждое «y» и каждый из соответствующих ( p + 1 ) -векторов «x» так, чтобы$p$$(p+1)$

$(y_{n+j}-x_{n+j}^T\beta)^2=\lambda\beta_j^2\,,\quad j=1,\ldots,p$

$y_{n+j}=0$$x_{n+j,j}=\sqrt{\lambda}$$x_{n+j,k}=0$$x_{n+j,0}=0$

потом

$(y_{n+j}-[x_{n+j,0}\beta_0+x_{n+j,1}\beta_1+x_{n+j,2}\beta_2+...+x_{n+j,p}\beta_p])^2=\lambda\beta_j^2$

Это работает для линейной регрессии. Это не работает для логистической регрессии, потому что обычная логистическая регрессия не минимизирует сумму квадратов невязок.

[Регрессия хребта - не единственное, что можно сделать с помощью таких уловок псевдонаблюдения - они встречаются в ряде других контекстов]

Обобщение этого рецепта для GLM действительно не сложно, поскольку GLM обычно подходят с использованием итеративно переоцененных наименьших квадратов . Следовательно, в каждой итерации можно заменить обычный шаг взвешенных наименьших квадратов на шаг взвешенных наименьших квадратов, наложенный на ребро, чтобы получить GLM, наложенный на ребро. Фактически, в сочетании с адаптивными штрафами гребня этот рецепт используется, чтобы соответствовать штрафованным L0 GLM (или лучшему подмножеству, то есть GLM, где штрафуется общее количество ненулевых коэффициентов). Это было реализовано, например, в пакете l0ara , см. Этот документ и этот для деталей.

Стоит также отметить, что самый быстрый способ решения правильной регрессии гребня в замкнутой форме использует

lmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  solve(crossprod(X) + diag(lambdas), crossprod(X, y))[, 1]
}

для случая, где n>=pили используя

lmridge_solve_largep = function (X, Y, lambda) (t(X) %*% solve(tcrossprod(X)+lambda*diag(nrow(X)), Y))[,1]

когда p>nи для модели без перехвата.

Это быстрее, чем использование рецепта увеличения строки , т.е.

lmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  qr.solve(rbind(X, diag(sqrt(lambdas))), c(y, rep(0, ncol(X))))
}

Если вам понадобятся ограничения неотрицательности для ваших установленных коэффициентов, тогда вы можете просто сделать

library(nnls)

nnlmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x
}

что затем дает чуть более точный результат, чем

nnlmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=rbind(X,diag(sqrt(lambdas))), b=c(Y,rep(0,ncol(X))))$x 
}

(и, строго говоря, только решение nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x является правильным).

Я еще не выяснил, как можно оптимизировать кейс с ограничением неотрицательности p > n- дайте мне знать, если кто-нибудь узнает, как это сделать ... [ lmridge_nnls_largep = function (X, Y, lambda) t(X) %*% nnls(A=tcrossprod(X)+lambda*diag(nrow(X)), b=Y)$xне работает]