第四章--R编程实践.doc

(word完整版)第四章 R编程实践
第四章 R编程实践

自此,我们已经对R 的功能有个全面的了解。一些数据分析都是在R对话窗口（R Console)中进行的，但对于复杂的数据分析显然是不行方便的。下面将从统计语言和编程角度来说明,我们将会了解一些在R编程实践中采用的简单的概念。
一、 循环和向量化
相比下拉菜单式的程序（R下开发的统计包R Commander)，R的一个优势在于它可以把一系列连续的操作简单的程序化.这一点上和所有其他计算机编程语言是一致的，但R有一些特性使得非专业人士也可以很简单的编写程序.
和其他编程语言一样,R有一些和C语言类似的控制结构.
（一)控制结构
1。条件语句：条件语句常用于避免除零或负数等数学问题.它有两种形式：
1)if (条件） 表达式1 else 表达式2
2）ifelse(条件，yes，no）
例如：
if （x>=0） sqrt(x） else NA
ifelse （x>=0,sqrt(x）,NA）
2。循环：有三种形式：
1)for （变量 in 向量) 表达式
2）while（条件） 表达式
3）repeat 表达式
例如:
假定我们有一个向量x，对于向量x 中值为b的元素，我们把0赋给另外一个等长向量y的对应元素,否则赋1。我们首先创建一个与向量x等长的向量y：
y 〈— numeric(length(x）)
for (i in 1:length（x)） if （x[i］ == b） y[i] <- 0 else y[i] <- 1
几个指令可以放在一个大括弧里面：
for (i in 1：length（x）） {
{
y［i］ 〈- 0
…
｝
x <— rnorm(10,—5，0.1）
y <- rnorm(10，5，2）
X <- cbind(x，y) ＃ 矩阵X的列名保持”x” 和”y"
for (i in 1：5/2) z［i] 〈— x[i］ + y[i］
另外一种可能的情况是当条件为真的时候一条指令才执行:
while (myfun > minimum） {
…
｝
count〈-0
while （count<10)｛
print（count)
count=count+1
}
#如果设置不正确while循环可以能导致无限循环，要小心
＃computeEstimate(）
x0〈-1
tol<—1e-5
repeat｛
x1<—rnorm(1,1,2)
if （abs(x1—x0)〈tol)｛
break
｝else{
x0〈-x1
}
}
#next,return
for（i in 1：100）{
if（i〈=20)｛
##skip
next
｝
＃do somthing here
｝
3.分支控制语句
switch（statement，list)
switch(2，’banana’，’apple’,'other’）
num〈—10
mode<-num%％3 ＃模3
cat('mode is',switch（mode+1,0，1,2),'\n’）
例子见“客观赋权法。R".
（二）向量化
但是，由于R的一些特性，在很多情况下循环和控制结构可以避免:向量化.向量化使得循环暗含在表达式中， 前面的例子中已经多次采用了.比如，两个向量之和：
> z 〈- x + y
这种加和可以通过循环结构来实现， 就像很多编程语言采用的策略一样:
> z <— numeric（length(x))
〉 for (i in 1:length(z)) z［i］ <- x[i］ + y[i]
在这个例子中,由于要用到向量的下标系统，有必要预先创建一个向量z.
显然，这种显式的循环仅仅用于向量x和y等长的情况:如果不是这样程序代码需要改变，而第一种表达方式在任何情况下都成立（如果对R的向量循环使用方式不了解的话，x和y不等长时要小心一点）。
条件语句(if… else） 可以用逻辑索引向量代替；同样上面的例子：
〉 y［x== b］ 〈— 0
〉 y[x!= b] <— 1
除了让代码更简洁，向量化的表达式在计算上效率更高，特别是大数据量的数据集.
（三)运用apply等函数
有多个’apply' 形式的函数用于避免使用代码的显式循环结构。apply 作用于矩阵的行或者/和列，它的语法规则是apply（X，MARGIN,FUN,…)，其中X是一个矩阵，MARGIN 表明是对行（1)还是列（2），或者二者(c（1，2)）进行操作，FUN 是一个函数（或操作符,但是这种情况下操作符要在一个括弧里面指定)，…是函数FUN可能的参数。一个简单的例子如下。
> x <- rnorm(10,-5，0.1）
> y 〈- rnorm(10,5，2)
> X <- cbind(x，y） # 矩阵X的列名保持”x" 和"y”
〉 apply(X，2，mean)
lapply（） 可以用于一个列表对象:它的语法类似apply 并且返回一个列表对象。
> forms <- list(y~x，y ～ poly（x,2)）
〉 lapply（forms,lm)
sapply（） 是函数lapply(） 一个更为灵活的变种，它可以接受向量或者矩阵作为主要参数，并且返回形式更为友好的结果，常常是表格方式。还有sweep(）、tapply（）函数等。
例如:
>attach（mtcars）
>tapply(mpg，cy1，mean）
〉col.mean=apply（mtcars，2，mean）
〉sweep（mtcars，2，col.mean,—)
表示从数据框mtcars中按列计算col.mean,并从mtcars中减去。
二、 用R写程序
一般情况下，一个R程序以ASCII 格式保存，扩展名为`.R’。如果一个工作要重复好多次，用R程序是一个不错的选择。在我们的第一个例子中，我们想对三个不同种属的鸟绘制一样的图，而且数据在三个不同的文件中.我们将一步一步的演示看R用不同的方式去完成这个简单的问题.
首先，我们凭直觉连续键入一系列命令，而且预先分割图形设备。
layout(matrix(1：3,3，1)) ＃分割图形界面
data 〈— read。table(”Swal.dat"） #读入数据
plot（data$V1，data$V2,type="l"）
title（”swallow"） ＃增加标题
data <- read。table（"Wren。dat”)
plot(data＄V1，data$V2，type=”l”）
title（"wren”)
data 〈- read。table(”Dunn。dat"）
plot（data＄V1，data＄V2，type="l”）
title("dunnock”)
字符`＃' 用于在程序中添加注释行：R会自动跳过注释行。
第一个程序的问题是在我们加入另外一个物种数据时,它过长。此外，一些命令多次执行， 因此它们可以放在一起，在执行的时候仅仅修改一些参数。这里的策略是把参数放到一个字符型的向量中去，然后用下标去访问这些不同的值。
layout(matrix（1：3,3，1)) ＃ 分割图形界面
species <— c（”swallow”,"wren","dunnock”)
file <— c（”Swal。dat"，"Wren。dat”,"Dunn。dat"）
for(i in 1：length（species）） {
data 〈- read。table(file［i］) ＃ 读入数据
plot（data＄V1，data＄V2，type=”l"）
title（species[i］) # 增加标题
}
注意代码read。table()里面的参数file［i] 上面没有双引号， 因为这个参数是字符型。
现在的代码比较紧凑.它比较容易加入新的物种， 因为设置物种名字和数据文件的向量都程序的前端.
如果扩展名为`。dat’的数据文件在R的工作目录下面，程序可以正常运行，否则用户要设置工作目录,或者设置绝对路径（例如：file <- ”/home/paradis/data/Swal。dat"）。如果程序保存在文件Mybirds.R 中，可以通过键入如下命令执行：
〉 source("Mybirds.R"）
和所有以文件作为输入对象的函数一样,如果该文件不在当前工作目录下面，用户需要提供该文件的绝对路径。
三、 编写你自己的函数
（一）编写自己的函数
1.函数的基本操作
大多数R的工作是通过函数来实现,而且这些函数的输入参数都放在一个括弧里面.用户可以编写他们自己的函数，并且这些函数和R里面的其他函数有一样的特性。
stdev=function(x)
{
＃计算标准差
sqrt（var(x)）
}
stdev
z=1：3
stdev（z)
函数的返回值可以是任何R对象,甚至可以是另一个函数。可以通过显式地调用return（）,把一个值返回给主调函数。如果不用这条语句，默认会将最后执行的语句的值作为返回值。例如:
oddcount<—function(x）{
k<—0#k的初始值
for(n in x){
if（n％%2==1)k=k+1#％％是模
｝
return(k)
｝
这个函数返回参数x中奇数的个数。
如果函数的输出包括多个变量,可以把他们集中在一个列表中，以列表为参数调用函数，让这个函数返回列表.
创建函数之后可以通过函数formals（）和body(）来查看其参数。
formals(oddcount）
body（oddcount）
可以使用page（）一个函数的代码,确定它的功能和用法.
page（abline）
函数内变量的赋值只在函数内有效。编写自己的函数可以让你有效的，灵活的与合理的使用R.我们再次使用前面读数据并且画图的例子。如果我们想在其他情况下进行这样的操作，写一个函数是一个不错的想法：
myfun <— function(S,F) ｛
data <— read。table（F）
plot(data$V1,data＄V2，type="l”）
title(S）
}
执行时,这个函数必须载入内存。当然，这有多种方式实现。和所有其他命令一样，函数的各行可以直接通过键盘键入，或者从一个文本编辑器里面拷贝粘贴。如果函数保存在一个文本文件中， 可以命令source(） 载入。如果用户期望一些函数在R启动时就被载入，可以把它们保存在工作目录下面的文件.RData中.另外一种方式是，配置文件`.Rpro¯le'或`Rpro¯le’ (详见？Startup for details)。最后,还可以创建一个包，但是这里不想多讨论.
一旦函数载入后，我们就可以键入一条命令以读入数据和画出我们想要的图，如myfun（”swallow","Swal.dat”）。因此，现在我们的程序有第三个实现的版本了：
layout（matrix(1:3,3，1))
myfun（"swallow","Swal。dat"）
myfun("wren”，"Wrenn。dat")
myfun（”dunnock”,"Dunn.dat”)
我们还可以用sapply（) 实现程序的第四个版本：
layout(matrix(1:3，3,1)）
species <- c("swallow"，"wren"，"dunnock”)
file <— c（"Swal.dat”，"Wren.dat”,”Dunn。dat"）
sapply（species，myfun,file)
2.环境与变量的作用域
在R里面，没有必要在一个函数里面进行变量声明。当一个函数执行时，R用一种称为词汇作用域（lexical scoping）的规则决定一个对象的作用域相对一个函数是局部还是全局。为了理解这种机制，我们可以认真研究一下下面的一个简单函数：
〉 foo <— function(） print(x）
> x <— 1
> foo（）
［1］ 1
x 不是为了在函数foo（）里面创建对象，因此R 将会在封装环境中搜索是否有个名为x 的对象， 和打印它的值（否则一条错误信息将会显示，执行中断）。
如果x是我们定义的函数中一个对象的名字,全局环境中变量x 值将会被采用。
> x <— 1
〉 foo2 〈— function(） { x <— 2；print(x) ｝
〉 foo2（)
[1］ 2
> x
［1］ 1
此时，print（) 使用在它所在的环境中定义的x,即环境foo2 中的x。
前面提及的”封装"是关键所在。在前面两个演示函数中，有两个环境:全局环境和函数foo 或foo2的局部环境。如果有三个或者更多的嵌套环境，对象搜索将逐层搜索直到全局环境。
rm（list=ls(））
ls(）
ls.str（)可以获得更多信息。
有两种方式指定一个函数的参数:通过它们的定义时的位置或者名字(又称为标签参数)。 例如，假定一个函数有三个参数：
foo 〈- function（arg1,arg2,arg3） {…}
foo（) 在执行时可以不用名字arg1,…,如果相应的参数对象放在相应的位置上,如：foo(x，y，z）。但是，如果使用了参数的名字，位置信息将会失效，如foo(arg3 = z，arg2 = y,arg1 = x）。R函数的另外一个特性是函数可能采用定义时的默认设置。例如：
foo <— function(arg1,arg2 = 5，arg3 = FALSE) ｛…}
命令foo（x）,foo(x，5，FALSE） 和foo（x，arg3 = FALSE) 将会得到一样的结果.使用一个函数的默认设置非常有用，特别在使用标签参数的时候，如foo(x,arg3 = TRUE） 仅仅改变一个默认设置.
环境层次中,某一层次的代码对它上级层次中所有变量至少有读的权限。但是,通过标准的运算符〈—直接对上级层次的变量进行写操作是不可行的。
#观察w的值是否有变化
w<—12
f<—function(y）{
d<—8
w<—w+1
y<—y-2
cat("w=",w，”\n”）
h〈-function(){
return(d＊（w+y）)
}
return(h（）)
｝
f(2)
cat(”w=”，w，"\n")
如果你希望对一个全局变量，或者更一般情况下，对当前层次的上级层次中任意变量进行写操作，可以在当前层次使用超赋值运算符〈〈-，或者函数assign（).
（二）函数的调用
函数也是对象,也可以给它们赋值，把它们用作其他函数的参数,如下所示：
＃函数的调用
f1<-function(a,b) return(a+b)
f2〈-function（a，b) return(a—b)
f<—f1
f(3,2)
f<—f2
f(3,2)
g<-function(h,a,b) h（a,b）
g(f1,3，2）
g(f2,3，2）
#在同一幅画上绘制若干个函数的图形
g1<—function(x) return（5*sin（x）+5）
g2<-function(x) return(sqrt(x^2+1）)
g3〈-function(x) return(2*x-1)
plot（c(0,6.5),c(—5，15)）#准备绘图，确定x和y的区域
for(f in c（g1，g2,g3））plot(f，0,6.5，add=T)#在原有的图形上绘图
（三)函数编辑
R提供了一个edit函数，让使用者可以在一个特定的窗口中更改自己的函数；fix使用预设的编辑器。
fix(stdev）
edit(file="for。r"）
我们来看另外一个例子.尽管这个例子不是纯的统计学例子，但是它很好地展示了R 语言的灵活性。假定我们想研究一个非线性模型的行为：这个模型(Ricker 模型）的定义如下：

这个模型广泛地用于种群动态变化（population dynamics）的研究里面，特别是鱼类的种群变化。 我们想用一个函数去模拟这个模型关于增长率r和初始群体大小N0的变化情况（承载能力K 常常设定为1且以这个值作为默认值）的影响；结果将以种群大小相对时间的图表示。我们还将设定一个可选项允许用户只显示最后若干步中种群大小(默认所有结果都会被绘制出来)。下面的函数就是做Ricker模型的数值模拟的。
ricker 〈- function(nzero，r,K=1,time=100,from=0,to=time） ｛
N 〈- numeric(time+1)
N［1］ 〈— nzero
for （i in 1:time） N［i+1] <— N［i]*exp(r＊(1 — N[i］/K）)
Time 〈- 0：time
plot(Time，N，type=”l"，xlim=c（from，to）)
}
你可以试一试下面的代码:
〉 layout（matrix（1:3,3,1)）
> ricker(0。01，0.1）;title（"r = 0。1”）
〉 ricker（0.05,0.2）；title（"r = 0。2")
〉 ricker（0。1，0。3)；title（”r = 0.3”）
四、 防错、除错与调试
1。防错与除错
R具备在函数中加入防错功能的应用函数：stop()和warning(）。stop()会传回错误信息并终止执行函数；warning()则在传回错误讯息的同时继续执行函数.
stdev=function（x）
{
＃计算标准差
if （！is。vector(x）)
stop("向量才能计算标准差”）
else sqrt（var（x)）
}
z=3
z=as。matrix（z)
stdev（z)
＃利用输出cat语句进行确认
g=function（x,y)｛x^2-y+1}
z〈—0
x〈—y^2+3＊g（z，2)
cat（"x=”,x，"\n"）
＃利用print语句进行调试
x<—y^2+3＊g（z,2）
cat(”x=”,x,"\n”）
w〈—28
q<—12
if(w+q>0）{
u<-1
print(”the ’if' was done”）
}else｛
v〈-10
print（"the ’else' was done”)
}
2。调试
R的核心调试工具由浏览器（browser）构成,它可以让你逐行运行代码，并在运行过程中检查。可以用有内置除错函数debug()或browser()函数来打开这一浏览器。
R的调试工具是针对单个函数的，如果你认为函数f（)中有漏洞，就可以调用debug(f）来设置函数的调试状态。调用undebug（f)将取消函数的调试状态。
debug(stdev）
stdev(z）
如果你只想对f（）进行一次调试会话，可以使用debugonce（f）函数，会在f（）第一次执行时对其设置调试状态。
你果你认为漏洞离函数的开始部分较远，可以在函数的某一行加入browser(）语句，那么浏览器只在程序执行到这一行时被打开。
#设置断点
f=function(x)｛x^3-2*x-1｝
root=function（x,y)｛
if(f(x)＊f(y）==0）browser（）
x1〈—x；x2〈-y
if (f（x1）＊f(x2）〉0)
stop（"区间端点的f（）值之正负号不相反”）
else{
if（f(x）>0） {x1〈-y;x2〈—x｝
｝
tol〈—1e-5
repeat｛
x0〈—(x1+x2）/2
if (f(x0）〉0){
x2<-x0
｝else{
x1〈-x0
｝
if （abs(x1-x2）〈tol）{
return（x1)
break
｝
}
｝
plot（f，xlim=c(-1.5，2))
abline（h=0)
root(-1,1)
进入调试浏览器时，命令提示符将从〉为变更为Browser［d]（d表示函数调用链的深度)。在命令提示符之后你可以输入以下命令：
n（代表next）：执行下一行然后暂停，直接键入Enter也一样.
c（代表continue)：与n类似,但是如果当前在某个循环中，将会执行剩余的循环退出循环才会暂停。如果你不在某个循环中,那么下一次暂停前函数的剩余部分将被执行。
任意的R命令：在浏览器中，你依然处在R的交互模式,因此可以任意查看变量的取值.
where：会输出一份栈跟踪路径，即显示到达当前位置的过程中函数的调用序列。
Q:该命令将会退出浏览器，返回R的交互模式.
setBreakpoint(“f.R”，10)

思考与练习题
4.1运行如下程序，观察其结果：
>total=0
〉for (i in c（1,2，4,5，9)）
>total=total+i
>total
4.2编写一个用二分法求非线性方程根的函数，并求方程：

在区间［1，2］内的根，精度要求e=10—5。
4.3根据“客观赋权法。R"中的程序，学习switch分支控制语句的用法(查找帮助文件）.
4。4。根据数据集hw1_data编写程序,运用apply、tapply等函数回答下面问题：
(1）数据集臭氧值（Ozone)、太阳光(Solar。R）、风力（Wind）和温度（Temp)的平均值各是多少？
（2）按月（Month）平均的 “Temp"值是多少？
4.5试编写函数is.sym(）检查矩阵是否为方阵,然后检查矩阵是否对称,输出True、False或适当的错误信息。
4.6试编写函数oddsevens()，返回一个参数x中奇数和偶数的个数。
4.7试编写函数findruns(）寻找由0、1构成的向量x连续出现长度为k的1的游程个数.例如对于向量（1,0，0,1，1，1，0,1，1），从它的第4个元素开始有长度为3的游程，而长度为2的游程分别开始于第4、5、8的位置。
4。8试编写函数pread（)进行天气预测(1代表有雨、0代表没有雨）.根据天气的历史数据x选择最近k天的天气记录预测下一天的天气；对历史数据进行训练(随机生成）,比较不同的k,选择误差最小的k用于未来的预测.
9