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

今天药忘吃喽~ 2022-09-16 08:50 222阅读 0赞

### CV岭参数的选择 ###

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

# 0引言 #

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

# 1、代码 #

## 1.1 获取数据函数 ##

**输入参数：**  
– p：变量的维数，一般取正整数，默认3维.（数值）  
– n: 大于p的一个整数，默认值为200.（数值）  
– Beta: 回归数据的真实参数，每个分量默认为1.（长度为p的向量）  
**输出参数：**  
– Data:  n ∗ ( p + 1 ) n\*(p+1) n∗(p\+1)的数据框，第一列为 y y y值，第二到 p + 1 p+1 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 ∗ ( p + 1 ) n\*(p+1) n∗(p\+1)的数据框，第一列为 y y y值，第二到 p + 1 p+1 p\+1列为p为自变量.  
– lambda: 岭回归的惩罚参数，默认值为0.1.（数值）  
**输出参数：**  
– beta：输出岭回归的参数估计.( p ∗ 1 p\*1 p∗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 ∗ ( p + 1 ) n\*(p+1) n∗(p\+1)的数据框，第一列为 y y y值，第二到 p + 1 p+1 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水平。水平有限发现错误还望及时评论区指正，您的意见和批评是我不断前进的动力。

[R]: https://blog.csdn.net/weixin_46111814/article/details/120564844?spm=1001.2014.3001.5501