绘制RSS

学习An Introduction to Statistical Learning with Applications in R (Second Edition) 时,第63页的FIGURE 3.2的图看起来挺有意思,如下:

 

 想试着自己画一下,花了不少时间终于画了个差不多

首先得到TV对应sales的回归函数

 1 > summary(ad)
 2        X                TV             radio          newspaper     
 3  Min.   :  1.00   Min.   :  0.70   Min.   : 0.000   Min.   :  0.30  
 4  1st Qu.: 50.75   1st Qu.: 74.38   1st Qu.: 9.975   1st Qu.: 12.75  
 5  Median :100.50   Median :149.75   Median :22.900   Median : 25.75  
 6  Mean   :100.50   Mean   :147.04   Mean   :23.264   Mean   : 30.55  
 7  3rd Qu.:150.25   3rd Qu.:218.82   3rd Qu.:36.525   3rd Qu.: 45.10  
 8  Max.   :200.00   Max.   :296.40   Max.   :49.600   Max.   :114.00  
 9      sales      
10  Min.   : 1.60  
11  1st Qu.:10.38  
12  Median :12.90  
13  Mean   :14.02  
14  3rd Qu.:17.40  
15  Max.   :27.00  
16 > tv_fit = lm(sales ~ TV, data=ad)
17 > tv_fit
18 
19 Call:
20 lm(formula = sales ~ TV, data = ad)
21 
22 Coefficients:
23 (Intercept)           TV  
24     7.03259      0.04754  

其中可以看到回归的model tv_fit的截距(b0)是7.03259, TV的斜率(b1)是0.04754

计算一下RSS

> deviance(tv_fit)
[1] 2102.531

> sum(resid(tv_fit)^2)
[1] 2102.531

用上面两种方式计算都可以,结果是一样的

下面就准备画图,首先我们准备好截距(b0)和斜率(b1)

> b0=seq(5,9,length=100)
> b1=seq(0.01,0.08,length=100)

参考FIGURE 3.2,b0取值范围5到9,b1取值范围0.01到0.08,各生成100个点

下面就有计算b0, b1对应的RSS的值了,为此定义了一个函数 rssloop,作为一个R语言的小白,不知道是否有更好的方法

> rssloop = function(b0,b1){
    ret = vector()
    len = length(b0)
    for(i in 1:len){
        ret = append(ret,sum((ad$TV*b1[i] + b0[i] - ad$sales)^2))
    }
    return(ret)
}

计算对应b0, b1范围的RSS

> f=outer(b0,b1,FUN=rssloop)

> dim(f)
[1] 100 100

f 是一个100 x 100的矩阵

 

然后调用contour函数画出等高线

> contour(b0,b1,f*0.001,levels=c(2.103,2.105,2.11,2.15,2.2,2.3,2.5,3), col="Blue",xlab="β0",ylab="β1")

画出的图显示如下:

 

 注意我把RSS的值除了个1000,FIGURE 3.2中应该也是除以1000后的值,但是原图中的红点不知道怎么画出来,不过可以看到趋势:等高线在趋近于一个点。

 使用下面的命令画出三维图

> persp(b0,b1,f,col="LightBlue")

 

 

感觉三维图和 FIGURE 3.2还是有些差别,奇怪,有时间再研究一下。

 

posted @ 2022-07-31 22:12  脱缰的野猪  阅读(46)  评论(0编辑  收藏  举报