Best subset ќе има најниско RSS на тренирачкото множество бидејќи со тоа ќе ги испробаме сите можни кобинации од k предиктори. Комбинациите добиени со другите два метода сигурно ќе бидат испробани и од best subset методот.
Не можеме да знаеме кој модел има најниско RSS на тестирачкото множество, бидејќи изборот на модел за k предиктори со кој било од овие методи се прави исклучиво земајќи предвид RSS на тренирачкото множество.
Точно, бидејќи forward stepwise ги зема предикторите од \(\mathcal{M}_k\) и избира еден од преостанатите \(p-k\) предиктори за да го добие \(\mathcal{M}_{k+1}\).
Точно, бидејќи backward stepwise ги зема предикторите од \(\mathcal{M}_{k+1}\) и избира еден од тие \(k+1\) предиктори кој треба да се отфрли за да го добие \(\mathcal{M}_k\).
Грешно, бидејќи двата метода комплетно независно еден од друг стигаат до резултатот и може да се многу различни.
Грешно, бидејќи двата метода комплетно независно еден од друг стигаат до резултатот и може да се многу различни.
Грешно, бидејќи изборот на \(k\) и \(k+1\) предиктори во best subset се комплетно независни еден од друг и може да се многу различни.
И со Lasso, и со Ridge, добиваме помалку флексибилен модел. Разликата е во тоа како се ограничени овие модели. Во двата случаи, моделот добива повисок bias, а му се намалува variance. Затоа, кога покачувањето на bias ќе биде пониско од намалувањето на variance, ќе добиеме подобра точност. Точниот одговор е под (iii).
Не линеарните методи се пофлексибилни и ја зголемуваат варијансата на моделот. Кога повишувањето на варијансата е пониско од паѓањето на bias, добиваме подобра точност. Точниот одговор е под (ii).
RSS на тренирачкото множество постепено ќе опаѓа, бидејќи моделот постепено ќе станува пофлексибилен. Точниот одговор е под (iv).
На почеток, RSS тестирачкото множество ќе опаѓа, бидејќи моделот ќе го намалува bias-от. Во одреден момент, моделот ќе стане премногу фелксибилен, и варијансата ќе расте повеќе од што се намалува bias-от, па RSS ќе почне да расте. Ќе се добие U форма. Точниот одговор е под (ii).
Како што расте s, варијансата постепено ќе се зголемува. Точниот одговор е под (iii).
Како што расте s, bias-от постепено ќе се намалува. Точниот одговор е под (iv).
Оваа грешка ќе е константна. Точниот одговор е под (v).
Како што ја зголемуваме ламбдата, се зголемува степенот на регуларизација, а со тоа се намалува флексибилноста (варијансата) на моделот. Следствено на тоа, RSS на тренирачкото множество постепено ќе расте. Точниот одговор е под (iii).
Очекувано е намалувањето на варијансата на почеток да го подобри моделот на тестирачкото множество (да го намали RSS), бидејќи ќе се намали over-fitting-от. После некоја вредност, RSS повторно ќе почне да расте, бидејќи намалувањето на варијансата ќе биде помало од наголемувањето на bias-от. Ќе се добие U форма. Точниот одговор е под (ii).
Варијансата постепено ќе опаѓа. Точниот одговор е под (iv).
Bias-от постепено ќе расте. Точниот одговор е под (iii).
Оваа грешка ќе е константна. Точниот одговор е под (v).
\[ \sum_{i=1}^{2}y_i=0\implies|y_1|=|y_2|\implies y_1^2=y_2^2=y^2\\ \sum_{i=1}^{2}x_{ij}=0\implies|x_{1j}|=|x_{2j}|\implies x_{1j}^2=x_{2j}^2=x_j^2\\ x_{i1}=x_{i2}=x_i\implies|x_{1}|=|x_{2}|\implies x_{1}^2=x_{2}^2=x^2\implies x_1+x_2=0\\ RSS= \sum_{i=1}^{2}(y_i-x_i\sum_{j=1}^{2}\hat\beta_j)^2=\\ \sum_{i=1}^{2}(y^2-2y_ix_i\sum_{j=1}^{2}\hat\beta_j+x^2\sum_{j=1}^{2}\sum_{t=1}^{2}\hat\beta_j\hat\beta_t)=\\ 2y^2-2(\hat\beta_1+\hat\beta_2)(y_1x_1+y_2x_2)+2x^2(\hat\beta_1+\hat\beta_2)^2\\ \]
\[ RSS\_RIDGE=RSS+\lambda(\hat\beta_1^2+\hat\beta_2^2)=\\ 2y^2-2(\hat\beta_1+\hat\beta_2)(y_1x_1+y_2x_2)+2x^2(\hat\beta_1+\hat\beta_2)^2+\lambda(\hat\beta_1^2+\hat\beta_2^2)\\ \begin{cases} \frac{\partial RSS\_RIDGE}{\partial \hat\beta_1}=-2(y_1x_1+y_2x_2)+4x^2(\hat\beta_1+\hat\beta_2)+2\lambda\hat\beta_1=0\\ \frac{\partial RSS\_RIDGE}{\partial \hat\beta_2}=-2(y_1x_1+y_2x_2)+4x^2(\hat\beta_1+\hat\beta_2)+2\lambda\hat\beta_2=0 \end{cases}\implies\\ \hat\beta_1-\hat\beta_2=0\implies\\ \hat\beta_1=\hat\beta_2 \]
Во двата случаи (\(sign(\hat\beta_1)=sign(\hat\beta_2)=1\) и \(sign(\hat\beta_1)=sign(\hat\beta_2)=-1\)) се добива систем на две идентични равенки, па има бесконечно многу решенија. Единствено е важно бетите да имаат ист знак.
\[ RSS\_LASSO=RSS+\lambda(|\hat\beta_1|+|\hat\beta_2|)=\\ 2y^2-2(\hat\beta_1+\hat\beta_2)(y_1x_1+y_2x_2)+2x^2(\hat\beta_1+\hat\beta_2)^2+\lambda(|\hat\beta_1|+|\hat\beta_2|)\\ \begin{cases} \frac{\partial RSS\_RIDGE}{\partial \hat\beta_1}=-2(y_1x_1+y_2x_2)+4x^2(\hat\beta_1+\hat\beta_2)+\lambda\frac{\hat\beta_1}{|\hat\beta_1|}=0\\ \frac{\partial RSS\_RIDGE}{\partial \hat\beta_2}=-2(y_1x_1+y_2x_2)+4x^2(\hat\beta_1+\hat\beta_2)+\lambda\frac{\hat\beta_2}{|\hat\beta_2|}=0 \end{cases}\implies\\ \frac{\hat\beta_1}{|\hat\beta_1|}-\frac{\hat\beta_2}{|\hat\beta_2|}=0\implies\\ sign(\hat\beta_1)=sign(\hat\beta_2) \]
Во двата случаи ја плотираме со црвено вертикалната линија што минува низ бетата која го дава минимумот според (6.13). Гледаме дека таа се сече со параболата во нејзиниот минимум.
y = 3
lambda = 1
beta = (-10000:10000) / 100
rss_ridge = (y - beta) ^ 2 + lambda * beta ^ 2
beta_best = y / (1 + lambda)
plot(beta, rss_ridge, type="l", col="green")
abline(v=beta_best, col="red")
y = 3
lambda = 1
beta = (-10000:10000) / 100
rss_ridge = (y - beta) ^ 2 + lambda * abs(beta)
beta_best = y - lambda / 2
plot(beta, rss_ridge, type="l", col="green")
abline(v=beta_best, col="red")
\[ Y_i=\beta_0+\sum_{j=1}^{p}\beta_iX_{ij}+err_i \implies err_i=Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij}\\ f(Y|X,\beta)=\prod_{i=1}^np(Y_i|X_i,\beta)=\prod_{i=1}^{n}p(err_i|Y_i,X_i,\beta)=\prod_{i=1}^n\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{err_i-0}{\sigma})^2}=\\ \prod_{i=1}^n\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij}}{\sigma})^2}=\\ (\frac{1}{\sigma\sqrt{2\pi}})^ne^{-\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2} \]
\[ p(\beta)=\frac{1}{2b}e^{\frac{-|\beta|}{b}}\\ p(\beta|X,Y)\propto f(Y|X,\beta)p(\beta)=(\frac{1}{\sigma\sqrt{2\pi}})^ne^{-\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2}\frac{1}{2b}e^{\frac{-|\beta|}{b}}=\\ =(\frac{1}{\sigma\sqrt{2\pi}})^n\frac{1}{2b}e^{-\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2-\frac{|\beta|}{b}} \]
Ова е исто како lasso со \(\lambda=\frac{2\sigma^2}{b}\).
\[ log(f(Y|X,\beta)p(\beta))=log((\frac{1}{\sigma\sqrt{2\pi}})^n\frac{1}{2b})-\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2-\frac{|\beta|}{b}\\ \underset{\beta}{\arg\max}[log(f(Y|X,\beta)p(\beta))]=\underset{\beta}{\arg\min}[\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2+\frac{1}{b}\sum_{j=1}^{p}|\beta_j|]=\\ \underset{\beta}{\arg\min}[\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2+\frac{2\sigma^2}{b}\sum_{j=1}^{p}|\beta_j|]=\\ \underset{\beta}{\arg\min}[RSS+\frac{2\sigma^2}{b}\sum_{j=1}^{p}|\beta_j|] \]
\[ p(\beta)=\frac{1}{c\sqrt{2\pi}}e^{-\frac{\beta^2}{2c}}\\ p(\beta|X,Y)\propto f(Y|X,\beta)p(\beta)=(\frac{1}{\sigma\sqrt{2\pi}})^ne^{-\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2}\frac{1}{\sqrt{2c\pi}}e^{-\frac{\beta^2}{2c}}=\\ =(\frac{1}{\sigma\sqrt{2\pi}})^n\frac{1}{\sqrt{2c\pi}}e^{-\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2-\frac{\beta^2}{2c}} \]
Ова е исто како ridge со \(\lambda=\frac{\sigma^2}{c}\).
\[ log(f(Y|X,\beta)p(\beta))=log((\frac{1}{\sigma\sqrt{2\pi}})^n\frac{1}{\sqrt{2c\pi}})-\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2-\frac{\beta^2}{2c}\\ \underset{\beta}{\arg\max}[log(f(Y|X,\beta)p(\beta))]=\underset{\beta}{\arg\min}[\frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2+\frac{1}{2c}\sum_{j=1}^{p}\beta_j^2]=\\ \underset{\beta}{\arg\min}[\sum_{i=1}^n(Y_i-\beta_0-\sum_{j=1}^{p}\beta_iX_{ij})^2+\frac{\sigma^2}{c^2}\sum_{j=1}^{p}\beta_j^2]=\\ \underset{\beta}{\arg\min}[RSS+\frac{\sigma^2}{c}\sum_{j=1}^{p}\beta_j^2] \]
set.seed(1337)
x = rnorm(100)
e = rnorm(100)
beta = c(3, 2, 4, 1)
y = beta[1] + beta[2] * x + beta[3] * x ^ 2 + beta[4] * x ^ 3 + e
plot(x, y)
library(leaps)
data = data.frame(x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10, y)
data
print_plot_best <- function(list, name, max=TRUE) {
print(name)
print(list)
if (max) {
best = which.max(list)
} else {
best = which.min(list)
}
print(paste("Best:", best))
plot(list, type="l", ylab=name)
points(best, list[best], col="red", pch=19)
}
get_best <- function(method) {
best_subset = regsubsets(y ∼ ., data, nvmax=10, method=method)
print(summary(best_subset))
print_plot_best(summary(best_subset)$rsq, paste("R2", method))
print_plot_best(summary(best_subset)$adjr2, paste("Adj. R2", method))
print_plot_best(summary(best_subset)$cp, paste("Cp", method), max=FALSE)
print_plot_best(summary(best_subset)$bic, paste("BIC", method), max=FALSE)
}
get_best("exhaustive")
Subset selection object
Call: get_best("exhaustive")
10 Variables (and intercept)
Forced in Forced out
x FALSE FALSE
x.2 FALSE FALSE
x.3 FALSE FALSE
x.4 FALSE FALSE
x.5 FALSE FALSE
x.6 FALSE FALSE
x.7 FALSE FALSE
x.8 FALSE FALSE
x.9 FALSE FALSE
x.10 FALSE FALSE
1 subsets of each size up to 10
Selection Algorithm: exhaustive
x x.2 x.3 x.4 x.5 x.6 x.7 x.8 x.9 x.10
1 ( 1 ) " " " " "*" " " " " " " " " " " " " " "
2 ( 1 ) " " "*" "*" " " " " " " " " " " " " " "
3 ( 1 ) "*" "*" "*" " " " " " " " " " " " " " "
4 ( 1 ) "*" "*" "*" "*" " " " " " " " " " " " "
5 ( 1 ) "*" "*" "*" " " "*" " " "*" " " " " " "
6 ( 1 ) "*" "*" "*" "*" " " "*" "*" " " " " " "
7 ( 1 ) "*" "*" "*" "*" "*" " " "*" "*" " " " "
8 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" "*" "*"
9 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" " "
10 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
[1] "R2 exhaustive"
[1] 0.8381714 0.9810520 0.9930514 0.9932037 0.9932949 0.9935209 0.9935242
[8] 0.9935423 0.9935445 0.9935446
[1] "Best: 10"
[1] "Adj. R2 exhaustive"
[1] 0.8365201 0.9806613 0.9928343 0.9929176 0.9929383 0.9931029 0.9930315
[8] 0.9929746 0.9928989 0.9928193
[1] "Best: 6"
[1] "Cp exhaustive"
[1] 2135.132722 167.236272 3.799762 3.700135 4.443067 3.327344
[7] 5.281478 7.032421 9.002628 11.000000
[1] "Best: 6"
[1] "BIC exhaustive"
[1] -172.9114 -382.7902 -478.5015 -476.1124 -472.8579 -471.6812 -467.1274
[8] -462.8016 -458.2299 -453.6277
[1] "Best: 3"
get_best("backward")
Subset selection object
Call: get_best("backward")
10 Variables (and intercept)
Forced in Forced out
x FALSE FALSE
x.2 FALSE FALSE
x.3 FALSE FALSE
x.4 FALSE FALSE
x.5 FALSE FALSE
x.6 FALSE FALSE
x.7 FALSE FALSE
x.8 FALSE FALSE
x.9 FALSE FALSE
x.10 FALSE FALSE
1 subsets of each size up to 10
Selection Algorithm: backward
x x.2 x.3 x.4 x.5 x.6 x.7 x.8 x.9 x.10
1 ( 1 ) " " " " "*" " " " " " " " " " " " " " "
2 ( 1 ) " " "*" "*" " " " " " " " " " " " " " "
3 ( 1 ) "*" "*" "*" " " " " " " " " " " " " " "
4 ( 1 ) "*" "*" "*" "*" " " " " " " " " " " " "
5 ( 1 ) "*" "*" "*" "*" " " "*" " " " " " " " "
6 ( 1 ) "*" "*" "*" "*" " " "*" "*" " " " " " "
7 ( 1 ) "*" "*" "*" "*" " " "*" "*" " " "*" " "
8 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " " "*" " "
9 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" " "
10 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
[1] "R2 backward"
[1] 0.8381714 0.9810520 0.9930514 0.9932037 0.9932385 0.9935209 0.9935213
[8] 0.9935323 0.9935445 0.9935446
[1] "Best: 10"
[1] "Adj. R2 backward"
[1] 0.8365201 0.9806613 0.9928343 0.9929176 0.9928788 0.9931029 0.9930284
[8] 0.9929638 0.9928989 0.9928193
[1] "Best: 6"
[1] "Cp backward"
[1] 2135.132722 167.236272 3.799762 3.700135 5.221208 3.327344
[7] 5.321841 7.169696 9.002628 11.000000
[1] "Best: 6"
[1] "BIC backward"
[1] -172.9114 -382.7902 -478.5015 -476.1124 -472.0197 -471.6812 -467.0822
[8] -462.6475 -458.2299 -453.6277
[1] "Best: 3"
get_best("forward")
Subset selection object
Call: get_best("forward")
10 Variables (and intercept)
Forced in Forced out
x FALSE FALSE
x.2 FALSE FALSE
x.3 FALSE FALSE
x.4 FALSE FALSE
x.5 FALSE FALSE
x.6 FALSE FALSE
x.7 FALSE FALSE
x.8 FALSE FALSE
x.9 FALSE FALSE
x.10 FALSE FALSE
1 subsets of each size up to 10
Selection Algorithm: forward
x x.2 x.3 x.4 x.5 x.6 x.7 x.8 x.9 x.10
1 ( 1 ) " " " " "*" " " " " " " " " " " " " " "
2 ( 1 ) " " "*" "*" " " " " " " " " " " " " " "
3 ( 1 ) "*" "*" "*" " " " " " " " " " " " " " "
4 ( 1 ) "*" "*" "*" "*" " " " " " " " " " " " "
5 ( 1 ) "*" "*" "*" "*" " " "*" " " " " " " " "
6 ( 1 ) "*" "*" "*" "*" " " "*" "*" " " " " " "
7 ( 1 ) "*" "*" "*" "*" " " "*" "*" " " " " "*"
8 ( 1 ) "*" "*" "*" "*" " " "*" "*" "*" " " "*"
9 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" " " "*"
10 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
[1] "R2 forward"
[1] 0.8381714 0.9810520 0.9930514 0.9932037 0.9932385 0.9935209 0.9935215
[8] 0.9935391 0.9935439 0.9935446
[1] "Best: 10"
[1] "Adj. R2 forward"
[1] 0.8365201 0.9806613 0.9928343 0.9929176 0.9928788 0.9931029 0.9930286
[8] 0.9929711 0.9928983 0.9928193
[1] "Best: 6"
[1] "Cp forward"
[1] 2135.132722 167.236272 3.799762 3.700135 5.221208 3.327344
[7] 5.319341 7.075975 9.010286 11.000000
[1] "Best: 6"
[1] "BIC forward"
[1] -172.9114 -382.7902 -478.5015 -476.1124 -472.0197 -471.6812 -467.0850
[8] -462.7527 -458.2213 -453.6277
[1] "Best: 3"
library(glmnet)
set.seed(1337)
x_data = data[, c(
"x", "x.2", "x.3", "x.4", "x.5", "x.6", "x.7", "x.8", "x.9", "x.10"
)]
y_data = data[, c("y")]
cross_validated = cv.glmnet(as.matrix(x_data), y_data, alpha=1)
plot(cross_validated)
print(paste("Best lambda:", cross_validated$lambda.min))
[1] "Best lambda: 0.114250899104086"
library(ISLR)
set.seed(1337)
train = sample(c(TRUE, FALSE), nrow(College), rep=TRUE)
test = (!train)
MSE на тренирачкото множество е 962395.8, MSE на тестирачкото множество е 1317253.
lin_model = lm(Apps ~ ., data=College[train, ])
train_preds = predict(lin_model, College[train, ])
test_preds = predict(lin_model, College[test, ])
sum((College[train, ]$Apps - train_preds) ^ 2) / nrow(College[train, ])
[1] 962395.8
sum((College[test, ]$Apps - test_preds) ^ 2) / nrow(College[test, ])
[1] 1317253
Со ламбда 301.6011, тест MSE изнесува 2166183.
library(leaps)
library(glmnet)
set.seed(1337)
grid = 10 ^ seq(10, -2, length=100)
X = model.matrix(Apps ~ ., data=College[train, ])
cross_validated = cv.glmnet(X, College[train, c("Apps")], alpha=0)
plot(cross_validated)
cross_validated$lambda.min
[1] 301.6011
ridge_reg = glmnet(
X,
College[train, c("Apps")],
alpha=0,
lambda=cross_validated$lambda.min
)
test_preds = predict(
ridge_reg,
model.matrix(Apps ~ ., data=College[test, ])
)
sum((College[test, ]$Apps - test_preds) ^ 2) / nrow(College[test, ])
[1] 2166183
Со ламбда 1.466566, тест MSE изнесува 1320929. Нема нули коефициенти.
set.seed(1337)
X = model.matrix(Apps ~ ., data=College[train, ])
cross_validated = cv.glmnet(X, College[train, c("Apps")], alpha=1)
plot(cross_validated)
cross_validated$lambda.min
[1] 1.466566
ridge_reg = glmnet(
X,
College[train, c("Apps")],
alpha=1,
lambda=cross_validated$lambda.min
)
test_preds = predict(
ridge_reg,
model.matrix(Apps ~ ., data=College[test, ])
)
sum((College[test, ]$Apps - test_preds) ^ 2) / nrow(College[test, ])
[1] 1320929
predict(
ridge_reg,
model.matrix(Apps ~ ., data=College[test, ]),
type="coefficients"
)
19 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) -8.748677e+02
(Intercept) .
PrivateYes -3.988067e+02
Accept 1.352464e+00
Enroll -8.724207e-01
Top10perc 4.831166e+01
Top25perc -1.492847e+01
F.Undergrad 1.651791e-01
P.Undergrad 4.791767e-02
Outstate -3.975135e-02
Room.Board 1.516512e-01
Books 1.725220e-01
Personal -4.353364e-03
PhD -5.240238e+00
Terminal -1.127289e+01
S.F.Ratio 2.638255e+01
perc.alumni -7.347307e+00
Expend 9.730295e-02
Grad.Rate 1.151520e+01
Користиме 17 компоненти. Тест MSE e 1317253.
library(pls)
set.seed(1337)
pcr_model = pcr(Apps ∼ ., data=College[train, ], scale=TRUE, validation="CV")
validationplot(pcr_model, val.type="MSEP")
test_preds = predict(pcr_model, College[test, ], ncomp=17)
sum((College[test, ]$Apps - test_preds) ^ 2) / nrow(College[test, ])
[1] 1317253
Користиме 17 компоненти. Тест MSE e 1317253.
set.seed(1337)
plsr_model = plsr(Apps ∼ ., data=College[train, ], scale=TRUE, validation="CV")
validationplot(plsr_model, val.type="MSEP")
test_preds = predict(pcr_model, College[test, ], ncomp=17)
sum((College[test, ]$Apps - test_preds) ^ 2) / nrow(College[test, ])
[1] 1317253
Нема некоја разлика во тест MSE, без разлика што користиме. Секогаш лошо се предвидуваат бројот на апликации.