和真实参数的均方误差 (MSE)之间的递减关系和估计值。$k$$\beta$$\hat{\beta}(k)$

该图显示了替代子集选择方法的结果。图片说明解释了实验设计：有 10 个非零其余 21 个元素为零。理想的子集选择方法将正确地报告哪些是非零的，哪些是零；换句话说，没有特征被错误地包含，也没有特征被错误地排除。$\beta$$\beta$$\beta$

当数据生成过程中的一个或多个特征被省略时，就会出现遗漏变量偏差。有偏参数估计的期望值不等于它们的真实值（这是偏差的定义），因此选择绘制说得通。（请注意，偏差的定义与这个实验设置并不完全一致，因为也是随机的。）换句话说，该图显示了对于各种子集选择方法当太小时（在这种情况下，当$\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$$\beta$$k$$k$$k<10$) 参数估计是有偏差的，这就是为什么图表显示对于小的值较大的原因。$\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$$k$

显然，情况不应该如此——向线性模型添加更多变量并不意味着更好地估计真实参数。

幸运的是，这不是情节所显示的。相反，该图显示，使用子集选择方法可以产生正确或不正确的结果，具体取决于的选择。$k$

然而，当添加额外的特征确实改善了参数估计时，这个图确实显示了一种特殊情况。如果建立一个表现出遗漏变量偏差的模型，那么包含这些变量的模型将实现较低的参数估计误差，因为不存在遗漏变量偏差。

添加更多变量确实意味着较低的训练误差，即较低的残差平方和。

您将本段中的演示与不使用子集选择的替代方案混淆了。一般来说，估计具有较大基数的回归会降低使用训练数据测量的残差；这不是这里发生的事情。

轴是否标记不正确？特别是，轴是否有可能显示残差平方和而不是？$y$$y$$\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$

我不这么认为；原始帖子中提出的推理路线本身并不能确定标签不正确。Sextus 的实验发现了类似的模式。它不完全相同，但曲线的形状足够相似。

顺便说一句，我认为由于该图显示 了实验的经验结果，因此根据 Cagdas Ozgenc 的建议，写出用于期望的估计量会更清楚。

ESL 中的图 3.6 是否正确？

回答这个问题的唯一确定方法是获取用于生成图形的代码。该代码未公开获得或由作者分发。

如果无法访问过程中使用的代码，则在标记图形或数据或系数的比例/位置时总是可能出现错误；Sextus 在使用标题中描述的程序重新创建图表时遇到问题的事实提供了一些间接证据，表明标题可能不完全准确。有人可能会争辩说，这些再现性问题支持标签本身或图形点可能不正确的假设。另一方面，描述可能不正确，但标签本身是正确的。

这本书的不同版本发布了不同的图像。但是不同图像的存在并不意味着任何一个都是正确的。

向线性模型添加更多变量并不意味着更好地估计真实参数

这不仅仅是估计变量，还有变量选择。当你只选择<10个变量时，你不可避免地会出错。

• 这就是为什么当您为子集选择较大尺寸时误差会减小。因为更多的系数，可能是来自真实模型的系数，正在被估计（而不是左等于零）。

• 由于变量之间的高度相关性，误差的减少比$k=10$

  最强的改进发生在 k=10 之前。但是在的情况下，您还没有到达那里，并且您偶尔会从真实模型中选择错误的系数。$k=10$

  此外，附加变量可能会产生一些正则化效果。

• 请注意，在某个点之后，大约 ，添加更多变量时错误会上升。$k=16$

图的再现

在最后的 R 代码中，我试图重现前向逐步情况的图表。（这也是这里的问题：Recreating figure from Elements of Statistical Learning）

我可以使图看起来相似

但是，我需要对生成进行一些调整，使用而不是（我仍然没有得到与开始于0.95 并下降到 0.65，而使用此处的代码计算的 MSE 则要低得多）。尽管如此，形状在质量上是相同的。$\beta \sim N(1,0.4)$$\beta \sim N(0,0.4)$

这张图中的误差并不是因为偏差：我想将均方误差分成偏差和方差（通过计算系数的平均误差和误差的方差）。但是，偏差非常低！这是由于参数之间的高度相关性。当您有一个只有 1 个参数的子集时，该子集中的选定参数将补偿缺失的参数（它可以这样做，因为它是高度相关的）。其他参数过低的量将或多或少地等于所选参数过高的量。因此，平均而言，一个参数或多或少会过高或过低。

• 上图的相关系数为 0.15 而不是 0.85。
• 此外，我使用了固定的和（否则偏差将平均为零，对此进行更多解释）。$X$$\beta$

参数估计的误差分布

下面你会看到参数估计中的误差是如何作为子集大小的函数分布的。这使得更容易理解为什么均方误差的变化表现得如此。$\hat\beta_1- \beta_1$

注意以下特点

• 小的子集大小有一个峰值。这是因为参数通常不包含在子集中，并且估计将为零，从而使误差等于。随着子集大小的增加和包含参数的概率增加，该峰值的大小减小。$\hat\beta$$\hat\beta - \beta$$-\beta$
• 当单峰尺寸减小时，或多或少有一个尺寸增大的高斯分布分量。这是参数包含在子集中时的错误。对于较小的子集大小，此组件中的误差不以零为中心。原因是该参数需要补偿另一个参数（与之高度相关）的遗漏。这使得偏差的计算实际上非常低。差异很大。

上面的例子是固定的和。如果您要更改每个模拟的，那么偏差每次都会不同。如果然后将偏差计算为，那么您将非常接近于零。$\beta$$X$$\beta$$\mathbb{E}(\hat \beta - \beta)$

library(MASS)

### function to do stepforward regression
### adding variables with best increase in RSS
stepforward <- function(Y,X, intercept) {
  kl <- length(X[1,])  ### number of columns
  inset <- c()
  outset <- 1:kl
  
  best_RSS <- sum(Y^2)
  ### outer loop increasing subset size
  for (k in 1:kl) {
    beststep_RSS <- best_RSS ### RSS to beat
    beststep_par <- 0
    ### inner looping trying all variables that can be added
    for (par in outset) {
      ### create a subset to test
      step_set <- c(inset,par)
      step_data <- data.frame(Y=Y,X=X[,step_set])
      ### perform model with subset
      if (intercept) {
        step_mod <- lm(Y ~ . + 1, data = step_data)
      }
      else {
        step_mod <- lm(Y ~ . + 0, data = step_data)
      }
      step_RSS <- sum(step_mod$residuals^2)
      ### compare if it is an improvement
      if (step_RSS <= beststep_RSS) {
        beststep_RSS <- step_RSS
        beststep_par <- par
      }
    }
    bestRSS <- beststep_RSS
    inset <- c(inset,beststep_par)
    outset[-which(outset == beststep_par)] 
  }
  return(inset)
}

get_error <- function(X = NULL, beta = NULL, intercept = 0) {
  ### 31 random X variables, standard normal 
  if (is.null(X)) {
    X <- mvrnorm(300,rep(0,31), M)
  }
  ### 10 random beta coefficients 21 zero coefficients
  if (is.null(beta)) {
    beta <- c(rnorm(10,1,0.4^0.5),rep(0,21))
  }
  ### Y with added noise
  Y <- (X %*% beta) + rnorm(300,0,6.25^0.5)
  
  
  ### get step order
  step_order <- stepforward(Y,X, intercept)

  ### error computation
  l <- 10
  error <- matrix(rep(0,31*31),31) ### this variable will store error for 31 submodel sizes
  for (l in 1:31) {
    
    ### subdata
    Z <- X[,step_order[1:l]]
    sub_data <- data.frame(Y=Y,Z=Z)
    
    ### compute model
    if (intercept) {
      sub_mod <- lm(Y ~ . + 1, data = sub_data)
    }
    else {
      sub_mod <- lm(Y ~ . + 0, data = sub_data)    
    }
    ### compute error in coefficients
    coef <- rep(0,31)
    if (intercept) {
      coef[step_order[1:l]] <- sub_mod$coefficients[-1]
    }
    else {
      coef[step_order[1:l]] <- sub_mod$coefficients[]
    }   
    error[l,] <- (coef - beta)
  }
  return(error)
}


### correlation matrix for X
M <- matrix(rep(0.15,31^2),31)
for (i in 1:31) {
  M[i,i] = 1
}

### perform 50 times the model 
set.seed(1)
X <- mvrnorm(300,rep(0,31), M)           
beta <- c(rnorm(10,1,0.4^0.5),rep(0,21)) 
nrep <- 500
me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### this line uses fixed X and beta
###me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### this line uses random X and fixed beta
###me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### random X and beta each replicate

### storage for error statistics per coefficient and per k
mean_error <- matrix(rep(0,31^2),31)
mean_MSE <- matrix(rep(0,31^2),31)
mean_var <- matrix(rep(0,31^2),31)

### compute error statistics
### MSE, and bias + variance for each coefficient seperately
### k relates to the subset size 
### i refers to the coefficient
### averaging is done over the multiple simulations
for (i in 1:31) {
  mean_error[i,] <- sapply(1:31, FUN = function(k) mean(me[k,i,]))
  mean_MSE[i,] <- sapply(1:31, FUN = function(k) mean(me[k,i,]^2))
  mean_var[i,] <- mean_MSE[i,] - mean_error[i,]^2
}


### plotting curves
### colMeans averages over the multiple coefficients
layout(matrix(1))
plot(1:31,colMeans(mean_MSE[1:31,]), ylim = c(0,0.4), xlim = c(1,31), type = "l", lwd = 2,
     xlab = "Subset size k", ylab = "mean square error of parameters",
     xaxs = "i", yaxs = "i")
points(1:31,colMeans(mean_MSE[1:31,]), pch = 21 , col = 1, bg = 0, cex = 0.7)
lines(1:31,colMeans(mean_var[1:31,]), lty = 2)
lines(1:31,colMeans(mean_error[1:31,]^2), lty = 3)

legend(31,0.4, c("MSE", "variance component", "bias component"),
       lty = c(1,2,3), lwd = c(2,1,1), pch = c(21,NA,NA), col = 1, pt.bg = 0, xjust = 1,
       cex = 0.7)

### plotting histogram
layout(matrix(1:5,5))
par(mar = c(4,4,2,1))
xpar = 1
for (col in c(1,4,7,10,13)) {
  hist(me[col,xpar,], breaks = seq(-7,7,0.05), 
       xlim = c(-1,1), ylim = c(0,500),
       xlab = "", ylab = "",         main=paste0("error in parameter ",xpar," for subset size ",col),
       )
}

这里有很好的答案，所以我会尽量保持简短并添加几点。

• 该图的重点是显示估计的斜率与其真实值的接近程度，而不是模型在样本外预测的程度，或者推断是否有效。$y$

向线性模型添加更多变量并不意味着更好地估计真实参数

• 不要认为这是添加更多变量。在所有情况下，您都是从先验确定的一组固定变量开始的。问题是您是否应该放弃其中一些变量来构建最终模型。根据您在数据中看到的内容删除变量通常是一件坏事。如果您保留所有变量（假设您有足够的数据，在这种情况下您这样做），您的估计将是无偏的。换句话说，在数据生成过程中斜率实际​​上为的变量在拟合模型中应该具有接近于它们应该大致正确。当您删除变量时，这不再一定是正确的。$0$$0$

  这种情况比较复杂，因为变量都是相互关联的。相关性意味着斜率与它们的真实值的变化比变量都相互正交时的变化更大。因此，如果您选择正确的变量，您可以在保持无偏性的同时稍微减少方差。然而...

我的直觉是 MSE 应该在最优$k$

• 那是因为您的直觉是逐步过程将选择正确的变量。不幸的是，这不一定会发生。您不太可能选择完全正确的变量。而且，如果您不只选择正确的变量，您将继续获得具有更高方差和有偏估计的抽样分布。

  现在，让我们考虑选择最好的，比如 15 或 20 个变量。我们将包含我们想要的 10 个并且只丢弃只会增加噪音的无价值变量的概率是多少？好多了。这就是为什么那里的曲线较低的原因。

因此，从中得出的结论是，如果您知道有多少变量是正确的，并且您知道它们都包含在您的数据集中，您可以专注于保留超出需要的部分比例，您可能只会扔掉垃圾. （当然，我觉得这些条件不太现实，而且这个讨论只涉及斜率估计，而不是样本预测或统计推断，所以我继续认为逐步程序是不明智的。）

它可以帮助您阅读网站上与这些主题相关的其他一些主题：

• 自动模型选择算法
• 为什么执行逐步选择后 p 值会产生误导？
• （您可能还想阅读我们分类在多重共线性标签。）

我尝试给出一个直观的答案，而没有实际检查并尝试重现代码。不知道图表是否错误，但我会解释它如何符合我的直觉。

问题是：“我认为它显示了子集大小 k 与真实参数 β 和估计值 β^(k) 的均方误差 (MSE) 之间的递减关系。显然，情况不应该如此 - 添加更多线性模型的变量并不意味着对真实参数的更好估计。（...）我的直觉是 MSE 应该在最佳 k 附近最低（由于相关性，在 5-10 之间）。

我认为正在发生的事情是这样的。这是关于变量选择的。如果选择了正确的 10 个变量，则估计 beta 的 MSE 应该是最小的。如果这些变量中的至少一个被遗漏，它应该会更大。请注意，相关性使这个问题变得更糟，因为如果遗漏了一个正确的非零 beta 变量，它的贡献将归因于由于相关性而已经在模型中的那些。这将使他们的估计器变得更糟，除了缺少的本身存在错误之外。它不是$\beta$确实，相关性的影响是，在估计量的 MSE 方面，我们可以用比 10 个正确的变量更少的变量做得很好。对于预测来说可能是正确的，因为缺失变量的信息被模型中已经存在的其他相关变量所补偿。但这不是图表的内容。可能对预测有帮助的相同效果将不利于估计，因为丢失的正确非零 beta 变量的影响将在模型中已经存在的变量中分配，从而影响它们的估计。

这意味着只有在总是或几乎总是准确地选择了 10 个正确变量时，最小值才会出现在 10 处。但这不太可能，因为相关性实际上使找到正确的变量变得非常困难。如果程序选择了 11、12 甚至 15 个变量，那么很可能会遗漏一个真正的非零 beta 变量。另一方面，真正的零贝塔变量可能无论如何都会具有相当低的估计系数，因此不会像错过正确的非零贝塔变量那样损害估计器 MSE。这在我看来解释了估计器 MSE 仅从甚至$k=16$$k=27$舞台上的左右。我觉得这一切都很好。它显示的是需要在此设置中选择多少变量才能以足够大的概率找到所有真正的非零值。16 对我来说似乎很现实，而且很明显 stagewise 在这个问题上遇到了困难，因为它需要许多步骤才能降低最初高估的参数。