绘制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还是有些差别,奇怪,有时间再研究一下。