统计模拟与R相关资料习题答案

1、统计模拟与,主讲： 朱 强 理学院,QQ群: 6052759 E-Mail:qiang_,Simulation and R,课程要求,1.本课程上机实习3-4次 2.每次上机需要完成一次实验并提交 3.本课程怎么考试？ 4.总评=考勤+实验作业+考试成绩,主要内容,1.R软件介绍 2.R统计分析 3.统计模拟介绍,教材和主要参考书籍,1、肖枝洪、朱强著, 统计模拟及其R实现， 武汉大学出版社，2010,2、Sheldon.M.Ross, Simulation,(3nd), Elsevier,2002,3、王兆军，陈广雷，邹长亮译，统计模拟(第4版), 人民邮电出版社，2007.,4、薛毅，陈立

2、萍，统计建模与R软件， 清华大学出版社,2006.,第一篇：R介绍,R是 一个开放(GPL)的统计编程环境 一种语言，是语言(由ATsexf 1 M F M M F Levels: F M 函数factor()用来把一个向量编码为一个因子，一般形式为 factor(x,levels=sort(unique(x),na.last=TRUE),labels, exclude=NA,ordered=FALSE) is.factor()检验对象是否是因子，as.factor()把向量转化为因子。levels()可以得到因子的水平 sex.level=levels(sexf);sex.level 1”F

3、” “M” 对因子向量可用table()统计各类数据的频数。 sex.tab=table(sexf);sex.tab,2.tapply()函数 例子：知道5位学生的性别，还知道这5位学生的身高，分组求身高的平均值 height=c(174,165,180,171,160) tapply(height,sex,mean) 1 F M 162.5 175.0 tapply() 的一般形式为 tapply(x,INDEX,FUN=NULL,simplify=TRUE) x是一对象，通常为一向量，INDEX是与X有同样长度的因子，FUN是要计算的函数 3.gl()函数 gl()可以方便地产生因子， g

4、l(n,k,length=n*k,labels=1:n,odered=FALSE) n为水平数，k为重复的次数，length为结果的长度，labels为n维向量，表示因子水平，ordered是逻辑变量，表示是否为有序因子，缺省值FALSE gl(2,3) gl(2,1,6),六、列表（list）和数据框（data.frame) 1、列表的构造 red=list(name=“黎明”,age=30,scores=c(85,76,90) rec $name 1 黎明 $age 1 30 $scores 1 85 76 90,列表元素可用”列表名下标”引用 rec2 130 rec32 但列表不同于向

5、量，每次只能引用一个元素，如rec1:2 是不允许的,对连续型变量产生因子,或分类 age- runif(100)*50 table(cut(age,c(0,10,20,30,40,50) (0,10 (10,20 (20,30 (30,40 (40,50 27 20 17 19 17,注意：”列表名下标”或”列表名下标范围”的用法是允许的，但与前面的含义不同，它还是一个列表 若指定了元素的名字，则引用列表元素还可以用它的名字作为下标 “列表名“元素名”或者“列表名$元素名” rec“age” rec$age 2、列表的修改 rec$name=“张三” rec$age=list(20,32,5

6、8) 可重新添加一些到列表中 rec$sex=“男” rec6=161 第5号元素没定义，所以值是”NULL”,若rec是一个向量，则空元素为”NA” 列表之间还可以连接起来， List.ABC=c(list.A,list.B,list.C),3.数据框 数据框通常是矩阵形式的数据，但矩阵各列可以是不同类型。数据框每列是一个变量，每行是一个观测。 (1)生成数据框 一般可以把数据框看作一种推广了的矩阵，它可用矩阵形式 表示，数据框可以用data.frame()函数生成，用法与list()函数相同 ls()#列出工作空间你定义的所有变量 str(x)#给出对象x的一些信息,d=data.fram

7、e(name=c(李明,张聪,王建),age=c(30, 35,28),height=c(180,162,175) d name age height 1 李明 30 180 2 张聪 35 162 3 王建 28 175 可以将列表转化为一个数据框:如比较他们之间的区别 d1=list(name=c(李明,张聪,王建),age=c(30, 35,28),height=c(180,162,175) as.data.frame(d1),（2）数据框的引用 引用数据框元素的方法与引用矩阵元素的方法相同。 d1:2,2:3 数据框的各变量也可用按列表引用(用双括号或$引用) d“age” 1 30

8、35 28 d$height 1 180 162 175 数据框的各行也可以定义名字，用rownames属性定义。 names(d) 1 name age height rownames(d)=c(one,two,three) rownames(d) 1 one two three,(3)attach()函数 数据框的主要用途是保存统计建模的数据。调用数据框的变量像以前d是不方便的,R提供了函数attach()可以把变量调入内存。 attach(d) r=age/height#此时age,height是可以利用的 但r的值不会进入数据框d，若想把刚才r当作新的数据框的变量 d$r=age/we

9、ight d name age height r one 李明 30 180 0.1666667 two 张聪 35 162 0.2160494 three 王建 28 175 0.1600000 detach(d)取消连接,（4）列表与数据框的编辑 可用函数edit()对列表或者数据框进行编辑 newd=edit(d) 当然也可对向量，数组，矩阵类型的数据进行修改和编辑. Fix 用函数merge(x, y, .) 对两个数据框合并成一个 可查看帮助.?merge,七、输入输出 1.输入 （1）read.table()等相关函数 read.table()函数是读表格形式的文件,若”班级”数据

10、已输入到一个纯文本文件”cjb.txt”中,假设存储在c:/cjb.txt, 其格式如下:,学号 姓名 平时 其中 期末 总评 07*001吴蒙蒙85666469 07*002高旭95809491 07*003李楠92785567 07*004赵利冲78484040 07*005何耀龙89585160 07*006何星85448174 07*007唐佳征95735365 07*008程西平96868990 07*009陈茂琴95846977 07*010程长久94798284 07*011王锋93748685 07*012何福顺90825366 07*013罗玲96939092 07*014刘

11、卿卿97679489 07*015张晓芳93728282 072301016杨乐遥88615060 072301017沈亮95415460 072301018高晶98769189 072301020郭慧93918387 072301021陈灿92766975 072301022李婷90766975 072301023陈元俊94745969 072301024谭海梅95646370 072301025张亚琴95898588 072301026胡阳98929897 072301027张冲92606569 072301028位会敏95928287 072301029胡红96789089 072301

12、030岳鑫璐96607979 072301031李欢欢988510097 072301032陈龙龙92656470 072301033朱盼88556668,读入格式为: cjb=read.table(c:/cjb.txt,header=TRUE) cjb is.data.frame(cjb) 1 TRUE header=TRUE表示第一行读入的是变量名,默认是 header=FALSE,即自动会添上列名.具体用法可见 帮助.?read.table d=read.csv(txt.csv)#读逗号分隔的文本文件 d=read.csv2(“txt.csv”)#读分号分隔的文本文件， 逗号代替小数点

13、d=read.delim(foo.txt)#读制表符分隔的文件 d=read.delim2(“foo.txt”)#读制表符分隔的文件， 逗号代替小数点 d=read.fwf(txt.fwf)#读固定宽的文件,D=readLines(“a.txt”)#每行有不同的结构 readLines(con = stdin(), n = -1L, ok = TRUE, warn = TRUE, encoding = unknown) cat(TITLE extra line, 2 3 5 7, , 11 13 17, file=ex.data, sep=n) readLines(ex.data, n=-1)

14、 unlink(“ex.data”) # 删除文件,R中一般不能直接读xls文件,可先转换为csv （2）scan()函数 scan()函数可直接读纯文本文件数据和直接从屏幕上 输数据 x=scan() 1: 2 4 4.2 5 5: Read 4 times x 1 2 4 4.2 5 假设有10名学生的体重数据已输入文件”weight.txt”, 格式为45 63.5 32.5 67.2 77.5 55 54 49 62 70 w=scan(c:/weight.txt),假设weight.txt中,奇数是女性体重,偶数是男性体重, 即数据中有不同的属性 wfm=scan(c:/weight

15、.txt,list(F=0,M=0) Read 5 records 将数据读入,并以列表的方式赋给变量wfm wfm $F 1 45.0 32.5 77.5 54.0 62.0 $M 1 63.5 67.2 55.0 49.0 70.0 其它数据格式的读入可以参见帮助文档中的R数据输入 输出.,2.输出 直接输入变量名即可输出 x=1:5 x 1 1 2 3 4 5 用print(x)与上式是等价的,里面有些参数,digits指定 每个数输出的有效数字位. options(digits=3) 1/3 1 0.333 write.table#输出数据框到文件 write.csv write.cs

16、v2 Write.matrix#输出矩阵到文件,需要加载包MASS,cat()函数也可用来输出,但它可以把多个参数连接起来 再输出(类似于paste()的功能) cat(i=,i,n)#注意为什么要用”n” cat(c(AB,C),c(E,F),n,sep=) ABCEF cat()还可以指定一个参数”file=文件名”,把结果写到 指定的文件中 cat(i=,1,n,file=c:/result.txt) 若指定的文件存在,则原内容被覆盖,但加上 append=TRUE可以不覆盖原文件而是在末尾附加.,八 程序控制结构,1.分支结构,if(条件) 表达式,或 if(条件)表达式1 else

17、表达式2 例如: if(any(x=0)y=log(1+x) else y=log(x),还有switch()函数,可参见讲义.,2.循环结构,(1)for循环,for(name in expr1)expr2 例如: 构造一个4阶Hilbert矩阵 n=4;x=matrix(0,nrow=4,ncol=4) for(i in 1:n) for(j in 1:n) xi,j=1/(i+j-1) , x ,1 ,2 ,3 ,4 1, 1.0000000 0.5000000 0.3333333 0.2500000 2, 0.5000000 0.3333333 0.2500000 0.2000000

18、3, 0.3333333 0.2500000 0.2000000 0.1666667 4, 0.2500000 0.2000000 0.1666667 0.1428571,(2)while循环,while(condtion)expr 当条件condtion成立,则执行表达式expr.例如编写一个计算100以内的Fibonacci数 f=1;f2=1;i=1 while(fi+fi+1f 11 1 2 3 5 8 13 21 34 55 89 144 13 233 377 610 987,(3)repeat循环,repeat expr repeat循环依赖break语句跳出循环.例如用repea

19、t循环编写一个计算100以内的Fibonacci数 f=1;f2=1;i=1 repeat fi+2=fi+fi+1 i=i+1 if(fi+fi+1=1000)break f 11 1 2 3 5 8 13 21 34 55 89 144 13 233 377 610 987,三个例子 myfn1- function(obs=10,n=100) x- rep(NA,100) for (i in 1:n) tmp - runif(obs) xi - mean(tmp) list(mn=mean(x),std=sd(x) Myfn1(),myfn2- function(obs=10) x - r

20、unif(obs) while(mean(x) 0.45) obs - 2*obs x - runif(obs) list(mn=mean(x),std=sd(x),obs=obs) myfn2(),myfn3= 0.45) break obs - 2*obs list(mn=mean(x),std=sd(x),obs=obs) myfn3(),九 编写函数,函数是一系列语句的组合，在R中可以写出自己的函数 形式为：变量名function(变量列表)函数体 函数的调用:变量名（变量值） 函数可以递归引用，但不提倡！ 例如： f=function(x) x+exp(x) f(3.5) 例子 使用

21、gamma函数求n!factorial = function(n) + if (n=0) gamma(n+1) + else print(“Please input a positive integer!) + factorial(6) factorial(-6) Choose(10,5) #求组合值,fzero=function(f,a,b,eps=1e-5) if(f(a)*f(b)0) list(fail=“finding root is fail!”) else repeat if(abs(b-a)eps)break x=(a+b)/2 if(f(a)*f(x)0)b=x else a

22、=x list(root=(a+b)/2,fun=f(x) ,例：编写一个用二分法求非线性方程根的函数，并求方程 在区间1,2内的根，精度要求,用刚才编写的函数计算上述方程的解 f=function(x)x3-x-1 fzero(f,1,2) $root 1 1.324718 $fun 1 -1.405875e-05,R软件提供了求一元方程的根的函数uniroot() uniroot(f,interval,lower=min(interval),upper=max(interval),tol=.Machie$double.eps0.25,maxiter=1000,.) uniroot(f,c(

23、1,2) $root 1 1.324718 $f.root 1 -5.634261e-07 $iter 1 7 $estim.prec 1 6.103516e-05,用于处理错误的函数 用于处理用户输入不正确的类型而可能出现的错误 warning() 若错误不严重以至影响整个计算 stop() 若错误可能导致计算中止 print() 显示必要的信息 formatC() 数值作为字符串输出 cat() 字符串联，可以插入n(换行)及t(tab键) paste() 字符粘贴(非字符型自动转换)例子：cat(R, is,a good, software.n) formatC(1/3, format

24、= f, digits = 4) formatC(1/3, format = e, digits = 4)paste(1:12) # 与as.character(1:12)等价paste(A, 1:6, sep = )paste(today is, date(),例.考虑某银行的存取系统. 设open.account()为存开户函数，deposit()为存款函数，withdraw()为取款函数， balance()为存取款平衡(余额)函数. 对于客户的存取过程如下：,open.account=function(total) list( deposit=function(amount) if(a

25、mount=0) stop(Deposits must be positive! n) total total) stop(You dont have that much money!n) total-total-amount cat(amount,withdrawn. Your balance is ,total,n n) , balance=function() cat( Your balance is ,total,nn) ) ,十、R在计算中的应用 1、求一元多项式方程的根 设一元多项式为 求方程f(x)=0的所有根。 例：求方程 的所有根。 a=c(-4,1,0,3) polyroo

26、t(a) 1 1.0-0.000000i -0.5+1.040833i -0.5-1.040833i,十、R在计算中的应用 2、求一元线性(非线性)方程的根 设一元函数为y=f(x) 求方程f(x)=0在区间a,b内的所有实根。 例：求方程 在区间(0,1)内的所有根。 两种方法：（1）先画出函数y=f(x)图像观察它的根的大致范围 f=function(x)(1-x)4-4*x*(1-x)3+4*x3*(1-x)-x4 x=seq(0,1,by=0.01);y=f(x) plot(x,y,type=l); abline(h=0,col=red),得到图像如下 看出在(0,1)间有三个根， 分

27、别位于(0.2,0.4), (0.4,0.6),(0.6，0.8) 利用函数uniroot求出其根 uniroot(f,c(0.2,0.4) uniroot(f,c(0.4,0.6) uniroot(f,c(0.6,0.8) 得到三个根分别为 0.211，0.5， 0.789 （2）或者直接利用包 rootSolve中的函数 uniroot.all()求出 (0,1)内的所有根 library(rootSolve) uniroot.all(f,c(0,1) 1 0.5000000 0.2113253 0.7886747,十、R在计算中的应用 3、求n元非线性方程组的根 例：求下列方程组的根 m

28、odel=function(x)c(F1=x12+x22-1,F2=x22-2*x1)#定义方程组 library(rootSolve) multiroot(model,c(1,1)#要先给初始值,$root 1 0.4142136 0.9101797 $f.root F1 F2 6.43565e-10 6.39019e-10 $iter 1 5 $estim.precis 1 6.41292e-10,十一、图形 R软件有很强的图形功能，可用简单的函数调用，迅速作出数据的各种图形。 比如：有一个班级的数据存贮在c:/work/class01.txt c1=read.table(c:/work/

29、class01.txt) colnames(c1)=c(Name,Sex,Age,Weight,Height) attach(c1)#建立连接，方便引用变量 plot(Height)画出身高的散点图 plot(Sex)画出性别的频数条形图,常用的绘图函数plot().也可作两个变量x与y的散点图 plot(Height,Weight,main=体重对身高的回归,xlab=身高,ylab=体重) plot(-50):50/25,(-50):50/25)3,type=l)表示画连线图 stem(Height)画茎叶图 boxplot(Weight)#绘制盒形图，可看出数据的大体分布，中间粗体线是中

30、位数，其它的有最小值，最大值，1/4分位数，3/4分位数,hist(Weight)#绘制直方图 fit1=lm(WeightHeight)#求体重对身高的回归 p1=predict(fit1,c1)对数据c1进行预测 boxplot(list(体重=Weight,预报=p1)画出盒形图 函数qqnorm(),qqline()提供了绘画正态QQ图和相应的直线的方法 qqnorm(Weight) qqnorm函数的用法 qqnorm(y,ylim,main=Normal Q-Q Plot，xlab=,Ylab=,plot.it=TRUE) qqline(Weight),Plot通常在一个窗口出现,

31、根据需要自动打开. 有时需要同时打开两个图形窗口,可用dev.new() 我们能够在窗口上保存图像,或者用如下命令保存 postscript(myfile.ps)#保存其他图形格#式的命令有jpeg(),bmp(),png(),tiff() plot(1:10) dev.off(),I.高水平绘图函数 1.plot() (1)plot(x,y)#x,y为向量，生成y关于x的散点图 (2)plot(x)#x为时间序列 (3)plot(f)#f为因子，生成f的直方图 plot(f,y)#f为因子,y为向量,生成y关于f水平的盒形图,(4)plot(df)#df为数据框 plot(expr)#exp

32、r为对象名称的表达式(如:a+b+c) plot(yexpr)#y为任意一个对象 例如: attach(c1) plot(c1)#三项指标构成的散点图 plot(Age+Height)#身高与年龄的散点图 plot(WeightAge+Height)#绘出两张散点图,一张是体重与年龄的,一张是体重与身高 (5)coplot()函数 coplot(WeightHeight|Age)#绘出按年龄段给出体重与身高的散点图,2.高水平绘图中的命令 (1)图中的逻辑命令 add=TRUE表示在所绘图在原图上加图,缺省为add=FALSE即新图替换旧图 axes=FALSE表示所绘图没有坐标轴,默认axe

33、s=TRUE (2)Type命令 type=“p”#绘制散点图(缺省值) type=“l”#绘制实线,type=“b” #画点和线 type=“o” #实线通过所有的点 type=“h” #绘出点到x轴的竖线 type=“s”or “S” #绘出阶梯形曲线 type=“n” #不绘任何点和曲线 (3)图中的字符串 axes=FALSE#不画坐标,默认 xlab=字符串,其字符串的内容是x轴的说明,ylab=字符串意义与前面相同,默认是变量名main=字符串,字符串内容是图的说明,出现在图的正上方. sub=字符串,副标题,出现在x轴的正下方,默认为空 xlim=c(lo,hi),ylim=c(

34、lo,hi)#坐标轴的范围,II.低水平作图函数 1.加点和线的函数 points()#在已有图形上加点,参数pch为画出指定的字符,可查看帮助.points(x,y)与plot(x,y)作用一样 lines()#在已有图上加直线,lines(x,y)与plot(x,y,type=“l”)类似,线的类型由参数lty=n决定,在help(par)中查看lty,2.在点处加标记 函数text()是在图上加标记,格式为 text(x,y,labels,)#x,y是数据形向量.labels可以是整数,也可是字符串,缺省是labels=1:length(x) 如: plot(Weight,Height,

35、type=n) text(Weight,Height) legend(x,y,legend,)#在制定位置给出一个盒子对图形进行解释,legend是标签字符串向量其他参数lty=,lwd=,col=,fill=,angle=,density,pch locator(n,type)#在鼠标点击的地方做标记,type 为做标记的类型(点或线),右键选终止 identify(x,y,labels)#在制定的点x,y处做标记,默认为1:length(x). 见demo4.2,3.在图上加直线 (1)abline(a,b)#表示画直线y=a+bx (2)abline(h=y)#表示画出一条过所有点的水平

36、直线 (3)abline(v=x)#表示画出一条过所有点的竖线 (4)abline(lm.obj)#表示画出线型模型得到的线型方程 (5)polygon(x,y,)#以数据(x,y)为坐标,依次连接所有的点, 绘出多边形,此时与plot(x,y,type=l)功能相似,4.在图上加标记,说明或其它内容 (1)title(main=“Main Title”,sub=“sub title”) 加上图的题目和子图的题目,主题目加在顶部,子题目加在底部 (2)axis(side,) side是边,side=1表示所加内容放在图的底部,side=2表示所加内容放在图的左侧,side=3表示所加内容放在图

37、的顶部 side=4表示所加内容放在图的右侧,5.其他类型的图 barplot()#画出一个条形图,可用一个向量指定条形的高度. boxplot()#产生一个箱型图 pie()#画饼图 pairs()#产生散点图的网格 contour(),persp(),image(),画2维图 见demo 4.3 在制定的地方可用数学标识 xlab=expression(hat(beta)2) 见demo4.4,6.一页多图,R可以在同一页面开若干个按行、列排列的窗格，在每个窗格中可以作一幅图。 每个图有自己的边空，而所有图的外面可以包一个“外边空”。 一页多图用mfrow参数或mfcol参数规定，如：,

38、par(mfrow=c(3,2) 表示同一页有三行两列共六个图，而且次序为按行填放。类似地， par(mfcol=c(3,2) 规定相同的窗格结构，但是次序为按列填放，即先填满第一列的三个再填第二列。 要取消一页多图只要再运行 par(mfrow=c(1,1)即可。 缺省时无外边空。为了规定外边空大小，可以用omi参数或oma参数。omi参数使用英寸为单位，oma参数以文本行高为单位，两个参数均为四个元素的向量，分别给出下、左、上、右方的边空大小。如： par(oma=c(2,0,3,0) 函数mtext用来在外边空加文字标注。其用法为 mtext(text, side = 3, line =

39、 0, outer = FALSE),其中text为要加的文本内容，side表示在哪一边写（1为下，2为左，3为上，4为右），line 表示边空从里向外数的第几行，最里面的一行是第0号，outer=TRUE时使用外边空，否则会使用当前图的边空。例如： par(mfrow=c(2,2), oma=c(0,0,3,0), mar=c(2,1,1,0.1) plot(Height) plot(Sex) boxplot(WeightSex);plot(Height,Weight) mtext(Simulation Data, outer=T, cex=1.5),在多图环境中还可以用mfg参数来直接跳到

40、某一个窗格，比如 par(mfg=c(2,2,3,2) 表示在三行两列的多图环境中直接跳到第二行第二列位置。mfg参数的后两个表示多图环境的行、列数，前两个表示要跳到的位置。 可以不使用多图环境而直接在页面中的任意位置产生一个窗格来绘图，参数为fig，如： par(fig=c(4,9,1,4)/10) 此参数为一个向量，分别给出窗格的左、右、下、上边缘的位置，取值为占全页面的比例，比如上面的例子在页面的右下方开一个窗格作图。,实验作业,1、从数据文件class.txt读入数据存为数据框，给数据框添加名字：name, age, height, weight, sex；在一页中对体重作出qqnor

41、m图、直方图、Boxplot； 对身高和体重作qqplot图，要求给出标题、坐标轴等。 2、对上述数据class.txt作尽可能多的本节讲过的图。要求同时用到高级图形和低级图形函数。,3、 将屏幕分割为四块，并分别画出 y=sin(x)， z=3*cos(x)， a=sin(x)*cos(x), b=sin(x)/x。,R-语言画图,x=0:18 plot(x,x,pch=x,col=x) points(x,18-x,pch=x) matplot(x,cbind(x,18-x),画图,spring= data.frame(compression=c(41,39,43,53,42,48,47,4

42、6),distance=c(120,114,132,157,122,144,137,141) attach(spring) (Hookes law: f=.5ks) plot(distance compression) plot(compression, distance),画图,par(mfrow=c(2,2)#准备画22的4个图 plot(compression, distance,main= Hookes Law) #只有标题的图 plot(compression, distance,main= Hookes Law, xlab= x,ylab= y) #标题+x,y标记 identif

43、y(compression,distance) #标出点号码,画图,plot(compression, distance,main=Hookes Law) #只有标题的图 text(46,120, f=1/2*k*s)#在指定位写入文字 plot(compression, distance,main=Hookes Law) #只有标题的图 text(locator(2), I am here!) #在点击的两个位置写入文字,画图,library(mass);data(Animals);attach(Animals) plot(body, brain) plot(sqrt(body), sqrt

44、(brain) plot(body)0.1, (brain)0.1) plot(log(body),log(brain),画图 (原始数据),par(mfrow=c(1,1), pch=1) plot(body, brain, xlim=c(0, 6000) text(x=body, y=brain,labels=s(Animals), adj=0)# adj=0 implies left adjusted text,画图(几个点),plot(bodyc(1,3,8), brainc(1,3,8),xlim=c(0,200) text(x=bodyc(1,3,8),y=brai

45、nc(1,3,8),labels=s(Animals,c(1,3,8), adj=0),画图 (对数变换后),Par(cex=0.7,mex=0.7) #character (cex) 语句之间用换行或分号(;)分开.,函数的编辑,fix(median.mean.ratio) 这时会出现一个记事本编辑器 编辑并试图存了之后如果发生错误(不会存),往往会让你编辑改过但有语法错误的fix(); 如果还使用fix(median.mean.ratio)则是编辑以前的版本.,函数,函数可以有缺省值,例如 Expo=function(y,x=2)z=yx;z y=seq(0,1,lengt

46、h=100);plot(y,Expo(y),type=l) 函数可以不写return,这时最后一个值为return的值.为了输出多个值最好使用list.,颜色画图函数,view.colours = function() plot(1, 1, xlim=c(0,14), ylim=c(0,3), type=n, axes=F, xlab=,ylab=) text(1:6, rep(2.5,6), paste(1:6), col=palette()1:6, cex=2.5) text(10, 2.5, Default palette, adj=0) rainchars = c(R,O,Y,G,B,

47、I,V) text(1:7, rep(1.5,7), rainchars, col=rainbow(7), cex=2.5) text(10, 1.5, rainbow(7), adj=0) cmtxt = substring(cm.colors, 1:9,1:9) # Split cm.colors into its 9 characters text(1:9, rep(0.5,9), cmtxt, col=cm.colors(9), cex=3) text(10, 0.5, cm.colors(9), adj=0) view.colours(),画图(继续),attach(spring);

48、spring s(spring)=LETTERS1:8 ;spring plot(distance,compression) identify(distance,compression,s(spring),画图(继续),par(mfrow=c(1,2) plot(x-0:50, y-pi*x2, xlab=Radius, ylab=expression(Area = pi*r2),type=b) plot(x-0:50, y-pi*x2,xlab=Radius, ylab=Area = pi*r2,type=l)#What is the difference?

49、par(mfrow=c(1,1),画图(继续),par(mfrow=c(2,2) plot(1:10,1:10,main=Straight line,type=l) locator(n=2,type=l)#加一条n个点决定的折线 hist(rnorm(50), main=Histogram of Normal) qqnorm(rt(100,2), main=Samples from t(2) plot(density(rnorm(50), main=Normal density) par(mfrow=c(1,1),画图(继续),plot(complex(arg=seq(-pi,pi,len=2

50、0),main=Complex numbers, sub=circle),画图(时间序列),x=runif(100);x=ts(x,start=1960+(3/12),frequency=12);options(digits=2);x plot(x,main=Time series),Jan Feb Mar Apr May Jun Jul Aug Sep 1960 0.28960 0.85129 0.05823 0.46995 0.57905 0.73799 1961 0.28672 0.20442 0.41045 0.96187 0.06345 0.18318 0.73739 0.08260

51、 0.74940 1962 0.43173 0.10568 0.44284 0.09334 0.44954 0.10872 0.05065 0.63924 0.31070 1963 0.67849 0.62165 0.21193 0.98964 0.96315 0.29877 0.95799 0.84184 0.69913 1964 0.99902 0.86650 0.83739 0.26195 0.70877 0.94352 0.56905 0.76416 0.95544 1965 0.15065 0.00066 0.87050 0.78968 0.58603 0.35572 0.04752

52、 0.25399 0.11059 1966 0.04642 0.74758 0.42476 0.50551 0.76536 0.73733 0.48250 0.07827 0.91807 1967 0.70602 0.94552 0.32180 0.20975 0.59760 0.04847 0.46288 0.40607 0.61738 1968 0.26741 0.01976 0.44851 0.90064 0.07154 0.59455 0.51154 Oct Nov Dec 1960 0.66690 0.06929 0.10811 1961 0.16218 0.59930 0.8723

53、5 1962 0.62558 0.57628 0.68994 1963 0.23256 0.93256 0.78502 1964 0.66011 0.79481 0.32449 1965 0.91167 0.44189 0.80639 1966 0.35219 0.89851 0.32350 1967 0.40915 0.51130 0.77735 1968,画图(matplot),sines=outer(1:20,1:4,function(x, y) sin(x / 20 * pi * y) matplot(sines, pch = 1:4, type = o, col = rainbow(

54、ncol(sines), options(digits=4) sines ,1 ,2 ,3 ,4 1, 1.564e-01 3.090e-01 4.540e-01 5.878e-01 2, 3.090e-01 5.878e-01 8.090e-01 9.511e-01 3, 4.540e-01 8.090e-01 9.877e-01 9.511e-01 4, 5.878e-01 9.511e-01 9.511e-01 5.878e-01 5, 7.071e-01 1.000e+00 7.071e-01 1.225e-16 6, 8.090e-01 9.511e-01 3.090e-01 -5.

55、878e-01 7, 8.910e-01 8.090e-01 -1.564e-01 -9.511e-01 8, 9.511e-01 5.878e-01 -5.878e-01 -9.511e-01 9, 9.877e-01 3.090e-01 -8.910e-01 -5.878e-01 10, 1.000e+00 1.225e-16 -1.000e+00 -2.449e-16 11, 9.877e-01 -3.090e-01 -8.910e-01 5.878e-01 12, 9.511e-01 -5.878e-01 -5.878e-01 9.511e-01 13, 8.910e-01 -8.09

56、0e-01 -1.564e-01 9.511e-01 14, 8.090e-01 -9.511e-01 3.090e-01 5.878e-01 15, 7.071e-01 -1.000e+00 7.071e-01 3.674e-16 16, 5.878e-01 -9.511e-01 9.511e-01 -5.878e-01 17, 4.540e-01 -8.090e-01 9.877e-01 -9.511e-01 18, 3.090e-01 -5.878e-01 8.090e-01 -9.511e-01 19, 1.564e-01 -3.090e-01 4.540e-01 -5.878e-01

57、 20, 1.225e-16 -2.449e-16 3.674e-16 -4.898e-16,画图(继续matplot),x - 0:50/50 matplot(x, outer(x, 1:8, function(x, k) sin(k*pi * x), ylim = c(-2,2), type = plobcsSh, main= matplot(,type = plobcsSh ),画图(符号和线条类型),par(mfrow=c(1,1) plot(1:10,sin(1:10),lty=1,pch=2,type=b),画图(legend),x - seq(-pi, pi, len = 65)

58、 plot(x, sin(x), type = l, ylim = c(-1.2, 1.8), col = 3, lty = 2) points(x, cos(x), pch = 3, col = 4) lines(x, tan(x), type = b, lty = 1, pch = 4, col = 6) title(legend(., lty = c(2, -1, 1), pch = c(-1,3,4), merge = TRUE), cex.main = 1.1) legend(-1, 1.9, c(sin, cos, tan), col = c(3,4,6), lty = c(2, -1, 1), pch = c(-1, 3, 4), merge = TRUE, bg=gray90),Table和barplot,tN - table(Ni - rpois(100, lambda=5);tN 0 1 2 3 4 5 6 7 8 9 11 12 1 7 12 14 16 1