第 1 章 常用函数及常见错误
1.1 产生数据
在进行模拟时, 我们经常会需要生成数据. 这里以正态分布为例, 说明如 何产生数据. 在 R 中, 每种分布都会有以下 4 个函数:
- 概率密度函数: dnorm(x, mean = 0, sd = 1, log = FALSE)
- 累计分布函数: pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
- 分位数函数: qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
- 随机数产生: rnorm(n, mean = 0, sd = 1):
其中前 3 个函数都支持向量输入, 即计算一组取值的概率密度、累计分布、 分位数. 最后一个函数常用来生成数据,n 即产生数据的个数. 如果需要生成 多维正态分布, 需要调用 MASS 包中的 *mvrnorm(n = 1, mu, Sigma, tol = 1e-6, empirical = FALSE, EISPACK = FALSE), 其中 Sigma 是指定的协方差矩阵.
1.2 定义运算符
空间模型 (SAR) 中, 会出现对角块矩阵 \[ \mathbf{W}=\left(\begin{array}{cccc}{\mathbf{M}} & {} & {} & {} \\ {} & {\mathbf{M}} & {} & {} \\ {} & {} & {\ddots} & {} \\ {} & {} & {} & {\mathbf{M}}\end{array}\right) \]
M 作为权重矩阵, 这种形式的矩阵可以利用克罗内克积 (Kronecker product, 符号为\(\otimes\)) 在 R 中很方便的产生,命令是%x%. 可以写成\[\mathbf{W} = \mathbf{I} \otimes \mathbf{M}.\]
diag(3)%x%matrix(1:6,2,3)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1 3 5 0 0 0 0 0 0
## [2,] 2 4 6 0 0 0 0 0 0
## [3,] 0 0 0 1 3 5 0 0 0
## [4,] 0 0 0 2 4 6 0 0 0
## [5,] 0 0 0 0 0 0 1 3 5
## [6,] 0 0 0 0 0 0 2 4 6这是借助自定义运算符实现的.自定义运算符是一种特殊的函数,当参数只有两个变量时,可以进行定义.用法如下:
'%myop%'<-function(a,b){a^b+b^a}
2%myop%3## [1] 17利用自定义运算符,可以实现很方便的功能.R中矩阵乘法(%*%)、Kronecker乘积(%x%)都是这样实现的.另外还有整除(%/%)和取余(%%)
9%/%4## [1] 213%%3## [1] 11.3 自动纠错
当我们输入的命令不规范时,R会自动纠正,以保证程序正常运行.
比如看下面的例子:
1:4 - 1:2## [1] 0 0 2 21:5-1:2## Warning in 1:5 - 1:2: longer object length is not a
## multiple of shorter object length## [1] 0 0 2 2 4当运算的向量长度不一致时,R会自动重复短的向量,使之长度与另外的向量长度相同进行运算.但是当长度是整数倍当时候,不会有任何提示.
再看下面当例子:
matrix(1:9,3,3)*(1:3)## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 4 10 16
## [3,] 9 18 27本来想计算矩阵与向量的乘积,但是错误使用了*, 未使用矩阵乘法%*%,R也可以计算返回一个矩阵,不会有warning.
当条件为向量是,只会判断第一个值:
if( 1 <= 1:3) print("真") else print("假")## Warning in if (1 <= 1:3) print("真") else print("假"):
## the condition has length > 1 and only the first element
## will be used## [1] "真"如果条件为向量,应该使用all或者any函数:
all(1 <= 1:3)## [1] TRUEall(2 <= 1:3)## [1] FALSE所有都真的时候返回TRUE.
any(4 <= 1:3)## [1] FALSEany(2 <= 1:3)## [1] TRUE只要有一个为真就返回TRUE.
1.4 predict函数
当我们拟合好一个模型时,下一步要做的就是评价模型好坏或者对新数据预测.这都需要将新的输入值带入模型中计算,得到预测值.区别只是有没有真值对比.对于简单模型,我们当然可以直接提取系数,自己计算预测值,但是当模型复杂时(比如时间序列模型),就不太容易操作.
R借助泛型函数1,编写模型都会提供summary、predict、plot等函数方便调用.但是在学习过程中发现经常会不小心错误使用,特此单独说明一下.下面以线性模型为例,先看正确的使用方法:
n=100;p=3;beta=c(1,2,3);
X = matrix(rnorm(n*p),n,p)
Y = X%*%beta + rnorm(n)
trainlist = sample(1:n,70)
regData = data.frame(Y,X)
fitmodel = lm(Y~.,data=regData[trainlist,])
pe = predict(fitmodel,newdata=regData[-trainlist,])
sum((pe-regData[-trainlist,1])^2)/length(pe)## [1] 0.9258关键点:
所有数据存在一个数据框中 通过下标控制训练集和测试集的数据
错误程序1:
n=100;p=3;beta=c(1,2,3);
X = matrix(rnorm(n*p),n,p)
Y = X%*%beta + rnorm(n)
fitmodel = lm(Y~X)
X2 = matrix(rnorm(n*p),n,p)
Y2 = X2%*%beta + rnorm(n)
#predict(fitmodel,newdata = X2)数据格式错误,不能执行
pe = predict(fitmodel,newdata = data.frame(X2))
sum((pe-Y2)^2)/length(Y2)## [1] 26.61这个程序能明显看出问题,误差不应该这么大,但是程序没有任何warning.
错误程序2:
n=100;p=3;beta=c(1,2,3);
X = matrix(rnorm(n*p),n,p)
Y = X%*%beta + rnorm(n)
fitmodel = lm(Y~X)
X2 = matrix(rnorm(0.2*n*p),0.2*n,p)
Y2 = X2%*%beta + rnorm(0.2*n)
pe = predict(fitmodel,newdata = data.frame(X2))## Warning: 'newdata' had 20 rows but variables found have
## 100 rowslength(pe)## [1] 100修改测试数据的条数,使之与训练集数据量不同,可以发现warning.提示我们数据行数不一样.并且我们测试集合X2是20行,但是预测值pe返回的是100个值.问题在于predict函数使用不正确.
我们用all指令查看:
all(predict(fitmodel)==predict(fitmodel,newdata = data.frame(X2)))## Warning: 'newdata' had 20 rows but variables found have
## 100 rows## [1] TRUE就是说我们输入的参数没起到作用.
查看函数说明2
?predict.lm| 参数 | 说明 |
|---|---|
| object | Object of class inheriting from “lm” |
| newdata | An optional data frame in which to look for variables with which to predict. If omitted, the fitted values are used. |
| … |
newdata不是必要的参数,当缺失时候会使用拟合模型的数据. 至于为什么输入的数据不能正确识别,因为名字不一样.R中参数的传递都是通过名字,我们在拟合lm时候,解释变量的名字叫‘X’,所以传入‘X2’不会识别到.如果把‘X2’改名成‘X’,即可正确预测3,比如看下面的程序:
n=100;p=3;beta=c(1,2,3);
X = matrix(rnorm(n*p),n,p)
Y = X%*%beta + rnorm(n)
fitmodel = lm(Y~X)
X = matrix(rnorm(0.2*n*p),0.2*n,p)
Y2 = X%*%beta + rnorm(0.2*n)
pe = predict(fitmodel,newdata = data.frame(X))
sum((Y2-pe)^2)/length(pe)## [1] 0.6031通过以上例子,想说明当程序的结果和预期不一致时.当然可能是我们的方法不对,但是也有可能是调用函数的方式出了问题.比如线性规划求解的lp函数,默认在正半轴求解4.
1.5 保存结果
我们重复了M次实验(运行几小时或者几天),计算了几个指标.但是后续通过阅读其他文献或者老师的建议,需要计算一个新的指标.这时候,就可以打开之前保存好的运行结果,只需要进行分析画图即可,不需要重新运行程序5.
getwd()
setwd()
save.image()getwd,setwd分别用于获取、设置当前的工作目录.有时候我们保存了结果,但是不知道存到哪里,可以通过getwd查看当前的工作目录.当然,最好是在程序执行前,手动设置好工作目录.
save.image用于保存工作空间,可以借助对于字符串操作函数paste等设置文件名.
1.6 结果输出
R只是我们模拟使用的工具,模拟结果需要以图表的形式在文章中展现.对于图片,只要单独保存成文件,在文章中插入即可.对于表格,如果不借助工具,会很耗时.这里我们通过xtable包中的函数,可以将表格数据转换成 LaTeX 内的表格形式.
这里顺便介绍一下R中package的安装和加载.不管是在R GUI中,还是在Rstudio中,都可以通过点选菜单进行安装.如果通过命令安装,如下
install.packages("xtable")安装后,并不能直接使用.需要通过library命令加载后才能使用.一个package中包含很多函数和数据,有时候我们只需要使用其中的一个函数,不需要加载整个包.这时候可以通过下面这样,直接调用某个包中的函数,xtable有一些参数可以设置输出的格式,比如下面指定了小数部分位数:
xtable::xtable(matrix(rnorm(15),3,5),digits=5)## % latex table generated in R 3.6.0 by xtable 1.8-4 package
## % Tue Aug 6 08:55:19 2019
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & 1 & 2 & 3 & 4 & 5 \\
## \hline
## 1 & -1.40045 & 0.78974 & -1.49898 & -0.49517 & 0.91428 \\
## 2 & 0.14569 & 0.76534 & 1.04113 & 0.68360 & 0.73087 \\
## 3 & -0.39239 & 0.15512 & -0.83092 & 0.77105 & -0.71431 \\
## \hline
## \end{tabular}
## \end{table}这只是最基础的表格,LaTeX中定制表头可以在这查看.