求助：完善R语言中岭回归准则计算函数的实现

Hey there! Let's get this ridge regression criterion function working properly. First, a quick recap: the standard ridge regression criterion combines the sum of squared residuals with a penalty term that shrinks the regression coefficients (usually excluding the intercept, since we don't penalize the constant term).

Issues with Your Current Code #

Let's break down what's missing or off in your initial draft:

• You're not initializing a variable to accumulate the sum of squared residuals (so the loop doesn't actually build up a total)
• The index beta[2:p+1] has a precedence problem in R—this evaluates to beta[(2:p)+1] instead of the intended beta[2:(p+1)]
• You're missing the critical ridge penalty term (λ * sum(coefficients²))
• Looping through each observation is inefficient in R; vectorized operations are faster and more idiomatic

Improved Vectorized Version (Recommended) #

This version uses R's vector/matrix operations to compute the criterion efficiently, which is way better for large datasets:

ridge.crit <- function(y, X, lambda, beta) {
  n <- length(y)
  p <- ncol(X)
  
  # Calculate residuals for all observations at once
  predicted <- beta[1] + X %*% beta[2:(p+1)]
  residuals <- y - predicted
  
  # Sum of squared residuals
  ss_residuals <- sum(residuals^2)
  
  # Ridge penalty (exclude intercept from penalty, standard practice)
  penalty <- lambda * sum(beta[2:(p+1)]^2)
  
  # Total criterion value
  final_answer <- ss_residuals + penalty
  
  return(final_answer)
}

If You Prefer a Loop (Not Recommended, But For Reference) #

If you want to keep the loop structure (for learning purposes), here's how to fix it:

ridge.crit_loop <- function(y, X, lambda, beta) {
  n <- length(y)
  p <- ncol(X)
  
  # Initialize sum of squared residuals to 0
  ss_residuals <- 0
  
  for(i in 1:n) {
    # Calculate predicted value for the i-th observation
    pred <- beta[1] + sum(beta[2:(p+1)] * X[i,])
    # Add the squared residual to our running total
    ss_residuals <- ss_residuals + (y[i] - pred)^2
  }
  
  # Add the ridge penalty term
  penalty <- lambda * sum(beta[2:(p+1)]^2)
  final_answer <- ss_residuals + penalty
  
  return(final_answer)
}

Test It Out! #

Let's verify with a sample dataset to make sure it works:

# Generate test data
set.seed(123)
n <- 100
p <- 3
X <- matrix(rnorm(n*p), nrow = n)
y <- 2 + 1.5*X[,1] - 0.8*X[,2] + 0.3*X[,3] + rnorm(n)
true_beta <- c(2, 1.5, -0.8, 0.3)
lambda <- 0.5

# Run the vectorized function
ridge.crit(y, X, lambda, true_beta)

This should return a value equal to the sum of squared noise terms (from the rnorm(n) we added) plus the penalty term 0.5*(1.5² + (-0.8)² + 0.3²).

内容的提问来源于stack exchange，提问作者John Thomas