glmnet 능형 회귀 분석이 왜 수동 계산과 다른 답변을 제공합니까? 몇 가지 결과를 얻었습니다.

능선 회귀 추정치를 계산하기 위해 glmnet을 사용하고 있습니다. 나는 glmnet이 실제로 내가 생각하는 것을하고 있다는 것을 의심하게 만드는 몇 가지 결과를 얻었습니다. 이를 확인하기 위해 solve에서 수행 한 능선 회귀와 glmnet의 결과를 비교하는 간단한 R 스크립트를 작성했습니다.

n    <- 1000
p.   <-  100
X.   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

beta1 <- solve(t(X)%*%X+5*diag(p),t(X)%*%Y)
beta2 <- glmnet(X,Y, alpha=0, lambda=10, intercept=FALSE, standardize=FALSE, 
                family="gaussian")$beta@x
beta1-beta2

차이의 규범은 일반적으로 20 정도이며 수치 적으로 다른 알고리즘으로 인해 발생할 수 없으므로 잘못된 일을해야합니다. glmnet릿지와 동일한 결과를 얻으려면 어떤 설정을해야 합니까?

답변

관찰하는 차이는 GLMNET이 목적 함수에서 사용하는 관측치 수 N의 추가 분할과 아래에 표시된 샘플 표준 편차에 의한 Y의 암시 적 표준화로 인한 것입니다.

\frac{1}{2 N} {‖ \frac{y}{s_{y}} - X β ‖}_{2}^{2} + λ ‖ β ‖_{2}^{2} / 2

$\frac{1}{2N}\left\|\frac{y}{s_y}-X\beta\right\|^2_{2}+\lambda\|\beta\|^2_{2}/2$

우리가 사용하는 곳 대신에, 에 대한 ,
$1 / n$

1 / n

$1/n$ $1 / (n - 1)$

1 / (n - 1)

$1/(n-1)$ $s_{y}$

s_{y}

$s_y$

s_{y} = \frac{\sum_{i} (y_{i} - \bar{y})^{2}}{n}

$s_y=\frac{\sum_i(y_i-\bar{y})^2}{n}$

베타에 대해 미분하고 방정식을 0으로 설정하면

X^{T} X β - \frac{X^{T} y}{s_{y}} + N λ β = 0

$X^TX\beta-\frac{X^Ty}{s_y}+N\lambda\beta =0$

베타를 해결하면 추정치를 얻습니다.

{\tilde{β}}_{G L M N E T} = (X^{T} X + N λ I_{p})^{- 1} \frac{X^{T} y}{s_{y}}

$\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}\frac{X^Ty}{s_y}$

Y의 원래 측정 항목에 대한 추정치 (및 해당 처벌)를 복구하기 위해 GLMNET은 추정치와 람다에 결과를 사용자에게 반환합니다. $s_{y}$

s_{y}

$s_y$

{\hat{β}}_{G L M N E T} = s_{y} {\tilde{β}}_{G L M N E T} = (X^{T} X + N λ I_{p})^{- 1} X^{T} y

$\hat{\beta}_{GLMNET}=s_y\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}X^Ty$

λ_{u n s t d .} = s_{y} λ

$\lambda_{unstd.}=s_y\lambda$

이 솔루션을 표준 능형 회귀 분석과 비교하십시오.

\hat{β} = (X^{T} X + λ I_{p})^{- 1} X^{T} y

$\hat{\beta}= (X^TX+\lambda I_p)^{-1}X^Ty$

것을 알 수 우리가 사용하는 경우, 또한 N의 추가 요인에 의해 조정됩니다 또는 기능, 페널티 암시 적으로 확장 할 것입니다 . 즉,이 함수를 사용하여 일부 대한 계수 추정치를 얻으면 대한 추정치를 효과적으로 얻는 입니다. $λ$

λ

$\lambda$ predict()coef() $1 / s_{y}$

1 / s_{y}

$1/s_y$ $λ^{*}$

λ^{*}

$\lambda^*$ $λ = λ^{*} / s_{y}$

λ = λ^{*} / s_{y}

$\lambda=\lambda^*/s_y$

이러한 관찰에 기초하여 GLMNET에 사용 된 페널티는 계수로 조정되어야합니다 . $s_{y} / N$

s_{y} / N

$s_y/N$

set.seed(123)

n    <- 1000
p   <-  100
X   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

beta1 <- solve(t(X)%*%X+10*diag(p),t(X)%*%(Y))[,1]

fit_glmnet <- glmnet(X,Y, alpha=0, standardize = F, intercept = FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

           [,1]        [,2]
[1,]  0.23793862  0.23793862
[2,]  1.81859695  1.81859695
[3,] -0.06000195 -0.06000195
[4,] -0.04958695 -0.04958695
[5,]  0.41870613  0.41870613
[6,]  1.30244151  1.30244151
[7,]  0.06566168  0.06566168
[8,]  0.44634038  0.44634038
[9,]  0.86477108  0.86477108
[10,] -2.47535340 -2.47535340

결과는 절편 및 표준화 된 X 변수의 포함으로 일반화됩니다. 표준화 된 X 행렬을 1의 열과 대각선 행렬을 포함하도록 수정하여 [1,1] 위치에 추가 영점을 갖습니다 (즉, 절편에 불이익을주지 않음). 그런 다음 각각의 표본 표준 편차를 기준으로 추정값을 표준화 해제 할 수 있습니다 (표준 편차를 계산할 때 1 / n을 사용하는지 확인).

{\hat{β}}_{j} = \frac{\tilde{β_{j}}}{s_{x_{j}}}

$\hat\beta_{j}=\frac{\tilde{\beta_j}}{s_{x_j}}$

{\hat{β}}_{0} = \tilde{β_{0}} - {\bar{x}}^{T} \hat{β}

$\hat\beta_{0}=\tilde{\beta_0}-\bar{x}^T\hat{\beta}$

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)
X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i]
}
X_scaled_ones <- cbind(rep(1,n), X_scaled)

beta3 <- solve(t(X_scaled_ones)%*%X_scaled_ones+1000*diag(x = c(0, rep(1,p))),t(X_scaled_ones)%*%(Y))[,1]
beta3 <- c(beta3[1] - crossprod(mean_x,beta3[-1]/sd_x), beta3[-1]/sd_x)

fit_glmnet2 <- glmnet(X,Y, alpha=0, thresh = 1e-20)
beta4 <- as.vector(coef(fit_glmnet2, s = sd_y*1000/n, exact = TRUE))

cbind(beta3[1:10], beta4[1:10])
             [,1]        [,2]
 [1,]  0.24534485  0.24534485
 [2,]  0.17661130  0.17661130
 [3,]  0.86993230  0.86993230
 [4,] -0.12449217 -0.12449217
 [5,] -0.06410361 -0.06410361
 [6,]  0.17568987  0.17568987
 [7,]  0.59773230  0.59773230
 [8,]  0.06594704  0.06594704
 [9,]  0.22860655  0.22860655
[10,]  0.33254206  0.33254206

인터셉트없이 표준화 된 X를 표시하는 코드가 추가되었습니다.

set.seed(123)

n <- 1000
p <-  100
X <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)

X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i]
}

beta1 <- solve(t(X_scaled)%*%X_scaled+10*diag(p),t(X_scaled)%*%(Y))[,1]

fit_glmnet <- glmnet(X_scaled,Y, alpha=0, standardize = F, intercept =
FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

             [,1]        [,2]
 [1,]  0.23560948  0.23560948
 [2,]  1.83469846  1.83469846
 [3,] -0.05827086 -0.05827086
 [4,] -0.04927314 -0.04927314
 [5,]  0.41871870  0.41871870
 [6,]  1.28969361  1.28969361
 [7,]  0.06552927  0.06552927
 [8,]  0.44576008  0.44576008
 [9,]  0.90156795  0.90156795
[10,] -2.43163420 -2.43163420

답변

https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html 에 따르면 가족이 gaussian인 glmnet()경우

\begin{matrix} (1) & \frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + λ \sum_{j = 1}^{p} (α | β_{j} | + (1 - α) β_{j}^{2} / 2) . \end{matrix}

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda\sum_{j=1}^p(\alpha|\beta_j| +(1-\alpha)\beta_j^2/2). \tag{1}$

사용시 glmnet(x, y, alpha=1)에 열 올가미 맞게 표준화를보고 된 위약금 용액 최소화하기위한 해결책이

그러나 적어도 능선 회귀 분석에 사용할 때보 고 된 페널티 에 대한 솔루션은 를 최소화하는 솔루션입니다

여기서 는 의 표준 편차입니다 . 여기서 페널티는 . $x$

x

$x$ $λ$

λ

$\lambda$

\frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + λ \sum_{j = 1}^{p} | β_{j} | .

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \sum_{j=1}^p |\beta_j|.$ glmnet_2.0-13glmnet(x, y, alpha=0) $λ$

λ

$\lambda$

\frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + λ \frac{1}{2 s_{y}} \sum_{j = 1}^{p} β_{j}^{2} .

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \frac{1}{2s_y} \sum_{j=1}^p \beta_j^2.$ $s_{y}$

s_{y}

$s_y$ $y$

y

$y$ $λ / s_{y}$

λ / s_{y}

$\lambda/s_y$

함수가 먼저 를 표준화 한 다음

는 효과적으로 을 최소화하는 것입니다.

또는 이와 동등하게
$y$

y

$y$ $y_{0}$

y_{0}

$y_0$

\begin{matrix} (2) & \frac{1}{2 n} \sum_{i = 1}^{n} (y_{0 i} - x_{i}^{T} γ)^{2} + η \sum_{j = 1}^{p} (α | γ_{j} | + (1 - α) γ_{j}^{2} / 2), \end{matrix}

$\frac{1}{2n} \sum_{i=1}^n (y_{0i}-x_i^T\gamma)^2 +\eta \sum_{j=1}^p(\alpha|\gamma_j| +(1-\alpha)\gamma_j^2/2), \tag{2}$

\frac{1}{2 n s_{y}^{2}} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + η \frac{α}{s_{y}} \sum_{j = 1}^{p} | β_{j} | + η \frac{1 - α}{2 s_{y}^{2}} \sum_{j = 1}^{p} β_{j}^{2},

$\frac{1}{2n s_y^2} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta \frac{\alpha}{s_y} \sum_{j=1}^p |\beta_j| +\eta \frac{1-\alpha}{2s_y^2} \sum_{j=1}^p \beta_j^2,$

\frac{1}{2 n} \sum_{i = 1}^{n} (y_{i} - β_{0} - x_{i}^{T} β)^{2} + η s_{y} α \sum_{j = 1}^{p} | β_{j} | + η (1 - α) \sum_{j = 1}^{p} β_{j}^{2} / 2.

$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta s_y \alpha \sum_{j=1}^p |\beta_j| +\eta (1-\alpha) \sum_{j=1}^p \beta_j^2/2.$

올가미 ( )의 경우 다시 스케일링 하여 가 의미가있는 것처럼 페널티를보고합니다 . 그런 다음 모든 에 대해 는 걸쳐 결과의 연속성을 유지하기위한 페널티로보고되어야 합니다. 이것은 아마도 위 문제의 원인 일 것입니다. 이것은 부분적으로 (2)를 사용하여 (1)을 해결하기 때문입니다. 또는 경우 에만 문제 (1)과 (2) 사이에 동등성이 있습니다 (즉, (1) 의 와 (2)의 간의 대응 ). 다른 모든 $α = 1$

α = 1

$\alpha=1$ $η$

η

$\eta$ $η s_{y}$

η s_{y}

$\eta s_y$ $α$

α

$\alpha$ $η s_{y}$

η s_{y}

$\eta s_y$ $α$

α

$\alpha$ $α = 0$

α = 0

$\alpha=0$ $α = 1$

α = 1

$\alpha=1$ $λ$

λ

$\lambda$ $η$

η

$\eta$ $α \in (0, 1)$

α \in (0, 1)

$\alpha\in(0,1)$ 문제 (1)과 (2)는 서로 다른 두 가지 최적화 문제이며 (1) 의 와 (2)의 사이에는 일대일 대응이 없습니다 . $λ$

λ

$\lambda$ $η$

η

$\eta$

How IT

언제든지 물어보세요.

glmnet 능형 회귀 분석이 왜 수동 계산과 다른 답변을 제공합니까? 몇 가지 결과를 얻었습니다.

답변

답변

답변