R语言入门——CV岭参数的选择

Posted 2021-11-06 统计学小王子

CV岭参数的选择

• 0引言
• 1、代码
• • 1.1 获取数据函数
  • 1.2 岭估计函数
  • 1.3 岭估计CV预测精度
  • 1.4 优化调用
• 2、总结

0引言

《R语言入门——多元回归交叉验证的实现》中介绍了使用交叉验证的思想使用RMSE衡量一个模型的预测精度，本文使用同样的思想确定岭回归的岭参数。具体的估计如下:
$$\hat{\beta }=(X^{T}X+\lambda E{)}^{-1}{X}^{T}Y.$$
本文就在于选择一个合适$\lambda$是的交叉验证后的RMSE最小，具体代码如下。

1、代码

1.1 获取数据函数

输入参数：
– p：变量的维数，一般取正整数，默认3维.（数值）
– n: 大于p的一个整数，默认值为200.（数值）
– Beta: 回归数据的真实参数，每个分量默认为1.（长度为p的向量）
输出参数：
– Data: $n\ast (p+1)$的数据框，第一列为$y$值，第二到$p+1$列为p为自变量.

getData <- function(n = 200, p = 3, Beta = rep(1, p)){
   library(MASS)
   x <- mvrnorm(n, rep(0, p), diag(rep(1, p)))
   names(x) <- paste0('x', 1:p)
   y <- x %*% Beta + rnorm(n, 0, 0.5)
   Data <- data.frame(y = y, x)
   Data
}

1.2 岭估计函数

输入参数：
– data: $n\ast (p+1)$的数据框，第一列为$y$值，第二到$p+1$列为p为自变量.
– lambda: 岭回归的惩罚参数，默认值为0.1.（数值）
输出参数：
– beta：输出岭回归的参数估计.($p\ast 1$的矩阵)

rlm <- function(data = getData(20,3), lambda = 0.1){
  library(MASS)
  x <- as.matrix(cbind(1, data[,-1]))
  y <- data[,1]
  n <- nrow(x); p = ncol(x)
  beta <- ginv(t(x) %*% x + lambda*diag(p)) %*% t(x) %*% y
  return(beta)
}

1.3 岭估计CV预测精度

输入参数：
– data: $n\ast (p+1)$的数据框，第一列为$y$值，第二到$p+1$列为p为自变量.（数据框）
– lambda: 岭回归的惩罚参数，默认值为0.1.（数值）
– k: 交叉验证的折数，默认值为5.（数值）
输出参数：（列表）
– yhat：输出岭回归的预测值.(列表中的第一个元素，为一个长度为n的向量)
– RMSE：输出岭回归的RMSE值.(列表中的第二个元素，为一个长度为1的向量)

CV.rlm <- function(data = getData(20,3), lambda = 0.1, k = 5, p = rep(1, k)/k){
  n <- nrow(data)
  e <- n - sum(round(n*p))
  ngroup <- round(n*p)
  if(e != 0){
    for(i in 1:e){
      ngroup[i] <- ngroup[i] + 1
    }
  } 
  start <- c(1, cumsum(ngroup)[1:(k-1)]+1)
  end <- cumsum(ngroup)
  ind <- sample(1:n)
  yhats <- c()
  for(i in 1:k){
    testdata <- data[ind[start[i]:end[i]],]
    traindata <- data[-ind[start[i]:end[i]],]
    testx <- as.matrix(cbind(1, testdata[, -1]))
    beta <- rlm(data = traindata, lambda = lambda)
    yhat <- testx %*% beta
    yhats <- c(yhats, yhat)
  }
  erro <- yhats - data$y[ind]
  RMSE <- sqrt((t(erro) %*% erro)/nrow(data))
  out <- list(yhat = yhats, RMSE = RMSE)
  return(out)
}

1.4 优化调用

> nlminb
function (start, objective, gradient = NULL, hessian = NULL, 
    ..., scale = 1, control = list(), lower = -Inf, upper = Inf) 

> CV.rlm(data = getData(20,3), lambda = 0.1)
$yhat
 [1]  2.5900983  2.7405474 -1.1402619 -0.5300552 -0.1751282  1.6769301
 [7]  0.3326216 -2.2282021  0.4628716  0.1187324  0.8532815 -0.5505310
[13]  1.1872715 -0.2164008 -0.2094051 -1.4429908  0.4309090  0.8978614
[19] -0.1471460  1.7677643

$RMSE
          [,1]
[1,] 0.6613792

> nlminb(1, function(x) CV.rlm(data = getData(20,3), lambda = x)$RMSE)
$par  # 结果
[1] 1.015625

$objective
[1] 0.4249641

$convergence
[1] 1

$iterations
[1] 2

$evaluations
function gradient 
      21        2 

$message
[1] "false convergence (8)"

2、总结

最后希望可以帮助大家提高R水平。水平有限发现错误还望及时评论区指正，您的意见和批评是我不断前进的动力。