R语言上机实验三
R语言综合实验报告
学号：2013310200629姓名：王丹学院：理学院专业：信息与计算科学成绩：日期：年月日基于工业机器人能否准确完成操作的时间序列分析摘要：时间序列分析是预测领域研究的重要工具之一,它描述历史数据随时间变化的规律,并用于预测数据。
本文首先介绍了一些常用的时间序列模型，包括建模前对数据的预处理、模型的识别以及模型的预测等。
通过多种方法分析所得到的数据，实现准确建模，可以得出正确的结论。
关键词:自回归(AR)模型，滑动平均(MA)模型，自回归滑动平均(ARMA)模型，ARMA最优子集一、问题提出，问题分析随着社会日新月异的发展，不断创新的科技为我们的生活带来了越来越多的便利。
机器人也逐渐走向了我们的生活，工厂里使用机器人去工作也可以大大减少生产成本，但为了保证产品质量，工厂使用的机器人应该多次测试，确保动作准确无误。
现有一批数据，包含了来自工业机器人的时间序列（机器人需要完成一系列的动作，与目标终点的距离以英寸为单位被记录下来，重复324次得到该时间序列），对于这些离散的数据，我们期望从中发掘一些信息，以便对机器人做更好的改进或者确定机器人是否可以投入使用。
但我们从中并不能看出什么，需要借助工具做一些处理，对数据进行分析。
时间序列分析是通过直观的数据比较或作图观测,去寻找序列中包含的变化规律,这种分析方法称为描述性时序分析。
在物理、天文、海洋学等科学领域,这种描述性时序分析方法经常能够使人们发现一些意想不到的规律,操作起来十分简单而且直观有效,因此从史前到现在一直被人们广泛使用,它也是我们进行统计时序分析的第一步。
我们将利用自回归(AR)模型、滑动平均(MA)模型以及自回归滑动平均(ARMA)模型去解决遇到的问题。
二、数据描述和初步分析下面是我们接收到的数据，数据来源：/~kchan/TSA.htm0.0011 0.0011 0.0024 0.0000 -0.0018 0.0055 0.0055 -0.00150.0047 -0.0001 0.0031 0.0031 0.0052 0.0034 0.0027 0.00410.0041 0.0034 0.0067 0.0028 0.0083 0.0083 0.0030 0.00320.0035 0.0041 0.0041 0.0053 0.0026 0.0074 0.0011 0.0011-0.0001 0.0008 0.0004 0.0000 0.0000 -0.0009 0.0038 0.00540.0002 0.0002 0.0036 -0.0004 0.0017 0.0000 0.0000 0.00470.0021 0.0080 0.0029 0.0029 0.0042 0.0052 0.0056 0.00550.0055 0.0010 0.0043 0.0006 0.0013 0.0013 0.0008 0.00230.0043 0.0013 0.0013 0.0045 0.0037 0.0015 0.0013 0.00130.0029 0.0039 -0.0018 0.0016 0.0016 -0.0003 0.0000 0.00090.0017 0.0017 0.0030 -0.0001 0.0070 -0.0008 -0.0008 0.00090.0025 0.0031 0.0002 0.0002 0.0022 0.0020 0.0003 0.00330.0033 0.0044 -0.0010 0.0048 0.0019 0.0019 0.0031 0.00200.0017 0.0014 0.0014 0.0039 0.0052 0.0020 0.0012 0.00120.0031 0.0022 0.0040 0.0038 0.0038 0.0007 0.0016 0.00240.0003 0.0003 0.0057 0.0006 0.0009 0.0040 0.0040 0.00350.0032 0.0068 0.0028 0.0028 0.0048 0.0035 0.0042 -0.0020-0.0020 0.0023 -0.0011 0.0062 -0.0021 -0.0021 0.0000 -0.0019-0.0005 0.0048 0.0048 0.0027 0.0009 -0.0002 0.0079 0.00790.0017 0.0034 0.0030 0.0025 0.0025 0.0004 0.0031 0.0057-0.0003 -0.0003 0.0006 0.0018 0.0022 0.0042 0.0042 0.0055-0.0005 -0.0053 0.0028 0.0028 0.0005 0.0036 0.0017 -0.0043-0.0043 0.0066 -0.0016 0.0055 -0.0011 -0.0011 -0.0049 0.00470.0056 0.0057 0.0057 -0.0002 0.0056 0.0037 0.0012 0.00120.0018 -0.0025 -0.0011 0.0027 0.0027 0.0039 0.0058 0.00030.0040 0.0040 0.0042 0.0000 0.0056 -0.0029 -0.0029 -0.00260.0016 0.0019 0.0015 0.0015 0.0007 0.0007 -0.0044 -0.0030-0.0030 0.0013 0.0029 -0.0010 0.0009 0.0009 -0.0016 0.00000.0000 0.0014 0.0014 -0.0003 0.0009 -0.0068 0.0003 0.0003-0.0012 0.0037 -0.0019 0.0023 0.0023 -0.0033 -0.0002 -0.00100.0021 0.0021 0.0026 -0.0002 0.0011 0.0028 0.0028 -0.00040.0026 -0.0015 0.0002 0.0002 0.0018 -0.0005 0.0004 -0.0008-0.0008 0.0018 0.0019 0.0029 -0.0022 -0.0022 0.0010 -0.00330.0020 0.0000 0.0000 0.0003 0.0007 -0.0009 -0.0035 -0.00350.0010 0.0007 0.0028 -0.0008 -0.0008 -0.0034 -0.0010 -0.0018-0.0021 -0.0021 -0.0006 -0.0018 -0.0046 -0.0017 -0.0017 -0.0001-0.0029 0.0020 -0.0049 -0.0049 -0.0021 -0.0027 -0.0018 -0.0015-0.0015 0.0051 -0.0002 0.0000 -0.0006 -0.0006 -0.0012 0.00120.0000 0.0021 0.0021 -0.0001 0.0022 0.0055 -0.0010 -0.00100.0048 0.0006 0.0026 0.0004 0.0004 0.0000 0.0000 0.00080.0044 0.0044 0.0002 0.0036这一群数目庞大的数据，以我们直观的判断，它们错综复杂，且毫无规律可言，根本不能从中得到有用的消息。
r语言编程实验报告总结
r语言编程实验报告总结
本次实验主要是对R语言编程的学习和掌握进行实践操作，通过实验了解R语言的基本语法和数据结构，掌握R语言的编程方法和数据分析技巧。
在实验中，我们学习了R语言的基础知识，如基本数据类型、变量、运算符、数据结构等。
同时，我们也学习了R语言的控制结构，如条件语句、循环语句等，这些控制结构可以帮助我们更好地控制程序的执行。
除此之外，我们还学习了R语言的函数和包的使用，在实验中我们使用了一些常用的包，如ggplot2包和dplyr包，这些包可以帮助我们更加方便地进行数据分析和绘图。
同时，我们也学习了如何自己编写函数，并且熟练掌握了函数的调用和参数传递。
通过实验，我们还学习了如何进行数据处理和数据分析，包括数据的读取和写入、数据的清洗和转换、数据的统计分析和可视化等等。
我们使用R语言对一些真实数据进行了处理和分析，这些数据包括房价、气温、人口等等。
在实验中，我们遇到了一些问题，如代码错误、数据异常等等，但是通过对问题的分析和解决，我们不断提升了自己的编程能力和数据分析技能。
综上所述，通过本次实验，我们深入了解了R语言的编程方法和数据分析技巧，掌握了一些常用的包和函数，并且在实践中熟悉了数据处理和分析的整个过程，这对我们今后的学习和工作都具有重要的
意义。
r语言实验报告
武夷学院实验报告课程名称：大数据挖掘与统计机器学习项目名称：多元统计分析与R语言建模姓名：__张树捷__专业：数学与应用数学班级：__2班_学号：__20171071203__同组成员：无_时间：______一、实验部分：P37 3m=seq(0,3000,by=300)hist(as.matrix(jie),m,freq=F,col=1:7)m=seq(0,3000,by=300)hist(as.matrix(jie),m,freq=T,col=1:7)m=seq(0,96000,by=3000)t<-as.matrix(jie)cumsum(t)hist(cumsum(t),m,freq=F,col=1:12,las=3) 累计频率图：p55 二、2barplot(apply(dier,1,mean),las=3)#按行作均值条形图barplot(apply(dier,2,mean))#按列作均值图条形图boxplot(dier)#按列作箱尾图boxplot(dier,horizontal = T)#箱尾图中图形按水平放置stars(dier,full=T,key.loc = c(13,1.5))#具有图例的360度星相图stars(dier,full=T,draw.segments = T,key.loc = c(13,1.5))#具有图例的3 60度彩色星相图summary(dier)#数据分析图表faces(dier,ncol.plot = 7)#按每行7个作脸谱图library(MSG)#加载msg包andrews_curve(dier)#绘制调和曲线图P87 4cor(disan)#计算相关系数矩阵plot(disan)#散点图library(psych)#加载psych包corr.test(disan)#y，x1，x2，x3的相关系数fm=lm(y~.,data=disan)#显示多元线性回归模型summary(fm)#多元线性回归系数t检验（拟合优度检验）由于P<0.05,于是在05.0=α水平处拒绝0H ，接受1H ，即x1,x2，x3与y 之间存在回归关系。
R语言实验三
实验三数组的运算、求解方程（组）和函数极值、数值积分【实验类型】验证性【实验学时】2 学时【实验目的】1、掌握向量的四则运算和内积运算、矩阵的行列式和逆等相关运算；2、掌握线性和非线性方程（组）的求解方法，函数极值的求解方法；3、了解 R 中数值积分的求解方法。
【实验内容】1、向量与矩阵的常见运算；2、求解线性和非线性方程（组）；3、求函数的极值，计算函数的积分。
【实验方法或步骤】第一部分、课件例题：1.向量的运算x<-c(-1,0,2)y<-c(3,8,2)v<-2*x+y+1vx*yx/yy^xexp(x)sqrt(y)x1<-c(100,200); x2<-1:6; x1+x22.x<-1:5y<-2*1:5x%*%ycrossprod(x,y)x%o%ytcrossprod(x,y)outer(x,y)3.矩阵的运算A<-matrix(1:9,nrow=3,byrow=T);AA+1 #A的每个元素都加上1B<-matrix(1:9,nrow=3); BC<-matrix(c(1,2,2,3,3,4,4,6,8),nrow=3); C D<-2*C+A/B; D #对应元素进行四则运算x<-1:9A+x #矩阵按列与向量相加E<-A%*%B; E #矩阵的乘法y<-1:3A%*%y #矩阵与向量相乘crossprod(A,B) #A的转置乘以Btcrossprod(A,B) #A乘以B的转置4.矩阵的运算A<-matrix(c(1:8,0),nrow=3);At(A) #转置det(A) #求矩阵行列式的值diag(A) #提取对角线上的元素A[lower.tri(A)==T]<-0;A #构造A对应的上三角矩阵qr.A<-qr(A);qr.A #将矩阵A分解成正交阵Q与上三角阵R的乘积，该结果为一列表Q<-qr.Q(qr.A);Q;R<-qr.R(qr.A);R #显示分解后对应的正交阵Q与上三角阵Rdet(Q);det(R);Q%*%R #A=Q*Rqr.X(qr.A) #显示分解前的矩阵5.解线性方程组A<-matrix(c(1:8,0),nrow=3,byrow=TRUE)b<-c(1,1,1)x<-solve(A,b); x #解线性方程组Ax=bB<-solve(A); B #求矩阵A的逆矩阵BA%*%B #结果为单位阵6.非线性方程求根f<-function(x) x^3-x-1 #建立函数uniroot(f,c(1,2)) #输出列表中f.root为近似解处的函数值，iter为迭代次数，estim.prec为精度的估计值uniroot(f,lower=1,upper=2) #与上述结果相同polyroot(c(-1,-1,0,1)) #专门用来求多项式的根，其中c(-1,-1,0,1)表示对应多项式从零次幂项到高次幂项的系数7.求解非线性方程组（1）自编函数： (Newtons.R)Newtons<-function (funs, x, ep=1e-5, it_max=100){index<-0; k<-1while (k<=it_max){ #it_max 表示最大迭代次数x1 <- x; obj <- funs(x);x <- x - solve(obj$J, obj$f); #Newton 法的迭代公式norm <- sqrt((x-x1) %*% (x-x1))if (norm<ep){ index<-1; break #index=1 表示求解成功}; k<-k+1 }obj <- funs(x);list(root=x, it=k, index=index, FunVal= obj$f)} # 输出列表（2）调用求解非线性方程组的自编函数funs<-function(x){ f<-c(x[1]^2+x[2]^2-5, (x[1]+1)*x[2]-(3*x[1]+1)) # 定义函数组J<-matrix(c(2*x[1], 2*x[2], x[2]-3, x[1]+1), nrow=2,byrow=T) # 函数组的 Jacobi 矩阵list(f=f, J=J)} # 返回值为列表 : 函数值 f 和 Jacobi 矩阵 Jsource("F:/wenjian_daima/Newtons.R") # 调用求解非线性方程组的自编函数Newtons(funs, x=c(0,1))8.一元函数极值f<-function(x) x^3-2*x-5 # 定义函数optimize(f,lower=0,upper=2) # 返回值 : 极小值点和目标函数f<-function(x,a) (x-a)^2 # 定义含有参数的函数optimize(f,interval=c(0,1),a=1/3) # 在函数中输入附加参数9.多元函数极值（1）obj <-function (x){ # 定义函数F<-c(10*(x[2]-x[1]^2),1-x[1]) # 视为向量sum (F^2) } # 向量对应分量平方后求和nlm(obj,c(-1.2,1))（2）fn<-function(x){ # 定义目标函数F<-c(10*(x[2]-x[1]^2), 1-x[1])t(F)%*%F } # 向量的内积gr <- function(x){ # 定义梯度函数F<-c(10*(x[2]-x[1]^2), 1-x[1])J<-matrix(c(-20*x[1],10,-1,0),2,2,byrow=T) #Jacobi 矩阵2*t(J)%*%F } # 梯度optim(c(-1.2,1), fn, gr, method="BFGS")最优点 (par) 、最优函数值 (value)10.梯形求积分公式（1）求积分程序： (trape.R)trape<-function(fun, a, b, tol=1e-6){ # 精度为 10 -6N <- 1; h <- b-a ; T <- h/2 * (fun(a) + fun(b)) # 梯形面积 repeat{h <- h/2; x<-a+(2*1:N-1)*h; I <-T/2 + h*sum(fun(x)) if(abs(I-T) < tol) break; N <- 2 * N; T = I }; I}（2）source("F:/wenjian_daima/trape.R") # 调用函数f<-function(x) exp(-x^2)trape(f,-1,1)（3）常用求积分函数f<-function(x)exp(-x^2) # 定义函数integrate(f,0,1)integrate(f,0,10)integrate(f,0,100)integrate(f,0,10000) # 当积分上限很大时，结果出现问题integrate(f,0,Inf) # 积分上限为无穷大ft<-function(t) exp(-(t/(1-t))^2)/(1-t)^2 # 对上述积分的被积函数 e 2 作变量代换 t=x/(1+x) 后的函数integrate(ft,0,1) # 与上述计算结果相同，且精度较高第二部分、教材例题：1.随机抽样（1）等可能的不放回的随机抽样:> sample(x, n) 其中x为要抽取的向量, n为样本容量（2）等可能的有放回的随机抽样:> sample(x, n, replace=TRUE)其中选项replace=TRUE表示有放回的, 此选项省略或replace=FALSE表示抽样是不放回的sample(c("H", "T"), 10, replace=T)sample(1:6, 10, replace=T)（3）不等可能的随机抽样:> sample(x, n, replace=TRUE, prob=y)其中选项prob=y用于指定x中元素出现的概率, 向量y与x等长度sample(c("成功", "失败"), 10, replace=T, prob=c(0.9,0.1))sample(c(1,0), 10, replace=T, prob=c(0.9,0.1))2.排列组合与概率的计算1/prod(52:49)1/choose(52,4)3.概率分布qnorm(0.025) #显著性水平为5%的正态分布的双侧临界值qnorm(0.975)1 - pchisq(3.84, 1) #计算假设检验的p值2*pt(-2.43, df = 13) #容量为14的双边t检验的p值4.limite.central( )的定义limite.central <- function (r=runif, distpar=c(0,1), m=.5,s=1/sqrt(12),n=c(1,3,10,30), N=1000) {for (i in n) {if (length(distpar)==2){x <- matrix(r(i*N, distpar[1],distpar[2]),nc=i)}else {x <- matrix(r(i*N, distpar), nc=i)}x <- (apply(x, 1, sum) - i*m )/(sqrt(i)*s)hist(x,col="light blue",probability=T,main=paste("n=",i), ylim=c(0,max(.4, density(x)$y)))lines(density(x), col="red", lwd=3)curve(dnorm(x), col="blue", lwd=3, lty=3, add=T)if( N>100 ) {rug(sample(x,100))}else {rug(x)}}}5.直方图x=runif(100,min=0,max=1)hist(x)6.二项分布B（10，0.1）op <- par(mfrow=c(2,2))limite.central(rbinom,distpar=c(10,0.1),m=1,s=0.9)par(op)7.泊松分布: pios（1）op <- par(mfrow=c(2,2))limite.central(rpois, distpar=1, m=1, s=1, n=c(3, 10, 30 ,50)) par(op)8.均匀分布：unif(0,1)op <- par(mfrow=c(2,2))limite.central( )par(op)9.指数分布：exp(1)op <- par(mfrow=c(2,2))limite.central(rexp, distpar=1, m=1, s=1)par(op)10.混合正态分布的渐近正态性mixn <- function (n, a=-1, b=1){rnorm(n, sample(c(a,b),n,replace=T))}limite.central(r=mixn, distpar=c(-3,3),m=0, s=sqrt(10), n=c(1,2,3,10)) par(op)11.混合正态分布的渐近正态性op <- par(mfrow=c(2,2))mixn <- function (n, a=-1, b=1){rnorm(n, sample(c(a,b),n,replace=T))}limite.central(r=mixn, distpar=c(-3,3),m=0,s=sqrt(10),n=c(1,2,3,10)) par(op)第三部分、课后习题：3.1a=sample(1:100,5)asum(a)3.2（1）抽到10、J、Q、K、A的事件记为A，概率为P(A)=1(5220)其中在R中计算得：> 1/choose(52,20)[1] 7.936846e-15（2）抽到的是同花顺P(B)=(41)(91) (525)在R中计算得：> (choose(4,1)*choose(9,1))/choose(52,5) [1] 1.385e-053.3#(1)x<-rnorm(1000,mean=100,sd=100)hist(x)#(2)y<-sample(x,500)hist(y)#(3)mean(x)mean(y)var(x)var(y)3.4x<-rnorm(1000,mean=0,sd=1) y=cumsum(x)plot(y,type = "l")plot(y,type = "p")3.5x<-rnorm(100,mean=0,sd=1) qnorm(.025)qnorm(.975)t.test(x)由R结果知：理论值为[-1.96,1.96]，实际值为：[-0.07929,0.33001]3.6op <- par(mfrow=c(2,2))limite.central(rbeta, distpar=c(0.5 ,0.5),n=c(30,200,500,1000))par(op)3.7N=seq(-4,4,length=1000)f<-function(x){dnorm(x)/sum(dnorm(x))}n=f(N)result=sample(n,replace=T,size = 1000)standdata=rnorm(1000)op<-par(mfrow=c(1,2)) #1行2列数组按列（mfcol）或行（mfrow）各自绘图hist(result,probability = T)lines(density(result),col="red",lwd=3)hist(standdata,probability = T)lines(density(standdata),col="red",lwd=3) par(op)。
统计模拟与R相关资料习题答案
第一篇：R介绍
– R是
• 一个开放(GPL)的统计编程环境 • 一种语言，是Ｓ语言(由AT&T Bell实验室的Rick Becker, John Chambers,Allan Wilks开发)的一 种方言(dialect) 之一，另一则为S-plus. • 一种软件，是集统计分析与图形直观显示于一体的统计 分析
• 参考资料 随软件所附pdf文档(help->manuals),随版 本更新: – W.N. Venables, D.M. Smith and the R DCT: Introduction to R -- Notes on R: A Programming Environment for Data Analysis and Graphics, 2003. /R web/Rnotes/R.html – R DCT, The R Environment for Statistical Computing and Graphics -- Reference Index,2003. – R DCT, R Data Import/Export, 2003. – R DCT, R Language Definition,2003 – R DCT, Writing R Extensions,2003
– R作为一个计划(project)，最早(1995年)是由 Auckland大学统计系的Robert Gentleman 和Ross Ihaka开始编制，目前由Ｒ核心开发小 组(R Development Core Team – 以后用R DCT表示)维护，他们完全自愿、工作努力负责， 并将全球优秀的统计应用软件打包提供给我们。我 们可以通过Ｒ计划的网站()了解有关Ｒ的最新信息和使用说明， 得到最新版本的Ｒ软件和基于Ｒ的应用统计软件包 .
实验 R语言..
patientID 11 22 33 44
age 25 34 28 52
diabetes status Type1 Poor Type2 Improved Type1 Excellent Type1 Poor
2 Basic data management 2.3 Data input
1. Entering data from the keyboard mydata <- data.frame() # Create an empty data frame (or matrix) mydata <- edit(mydata) # invoke the text editor
1 Introduction to R 1.2 Obtaining and installing R
The Comprehensive R Archive Network (CRAN) at：
1 Introduction to R 1.3 A simple R session
2 Basic data management 2.4 Data management
fix(mydata)
# load and save data
2 Basic data management 2.3 Data input
2. Importing data from a delimited text file grades <- read.table("studentgrades.csv", header=TRUE, sep=",") grades
1 Introduction to R 1.4 Getting help
help.start( ) # general help help("mean") or ?mean # help on the function
R语言回归分析和方差分析上机
1.one-way ANOVA is.factor(wghtcat) [1] TRUE > levels(wghtcat) [1] "< 140" "> 200" "140-170" "170-200" > is.factor(agec) [1] TRUE > levels(agec) [1] "35-40" "41-45" "46-50" "51-55" "56-60"
new2=data.frame(age=60) lm.pred2=predict(M2,new2,interval="conf",le vel=0.95) print(lm.pred2) fit lwr upr 134.8568 133.4627 136.251 new3=data.frame(weight=160,age=60) lm.pred3=predict(M3,new3,interval="conf",le vel=0.95) print(lm.pred3) fit lwr upr 133.3425 131.9832 134.7018
new1=data.frame(weight=160) lm.pred1=predict(M1,new1,interval="conf",le vel=0.95) print(lm.pred1) fit lwr upr 126.8264 126.2617 127.3911
by(sbp, wghtcat, FUN=meansd)
r软件 实训内容
r软件实训内容
R软件是一种用于统计计算和可视化的编程语言和软件环境。
以下是一些可能的R软件实训内容：
1. R语言基础：学习R语言的语法、变量、数据结构、控制流、函数等基础知识。
2. 数据导入和整理：学习如何从各种数据源导入数据到R中，并掌握数据清洗、数据转换和数据重塑等技术。
3. 数据可视化：学习如何使用R中的各种可视化包（如ggplot2、lattice 等）创建各种图表和图形，以更好地理解数据。
4. 统计分析：学习如何使用R进行各种统计分析，如描述性统计、参数检验、非参数检验、回归分析、方差分析等。
5. 机器学习和数据挖掘：学习如何使用R中的各种机器学习包（如caret、randomForest、e1071等）进行数据挖掘和机器学习。
6. 实践项目：通过实际项目，将所学的R知识和技能应用到实际问题中，提高解决实际问题的能力。
以上是一些常见的R软件实训内容，具体实训内容可能会根据课程要求和实际需求而有所不同。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
理学院实验报告
班级：学号：姓名：实验编号：
实验三：概率和分布的R实现
一、实验目的与要求：
１、会用R给出常见分布的概率密度、概率、分位数和随机数。
２、会利用sample命令进行随机抽样，prod，choose命令计算概率。
3、会利用R绘制各类分布的图形。
4、会利用choose，prod命令计算古典概率。
二、实验内容：
1.从一副扑克牌(52张)中随机抽5张，求下列概率
(1) 抽到的是10，J，Q，K，A；
> 4/choose(52,5)
[1] 1.539077e-06
(2) 抽到的是同花顺。
> 9*choose(4,1)/choose(52,5)
[1] 1.385169e-05
注：同花顺是指5张同一色牌能按从小到大连续排序，如2<3<4<5<6，3<4<5<6<7，…，10<J<Q<K<A都为同花顺。
2.模拟随机游动:
(1)从两点分布中产生1000个随机数;
> x<-rbinom(1000,1,0.5)
> x
(2)用函数ifelse( )将上面随机数中的0替换成-1;
> ifelse(x==0,-1,1 )
(3)用函数cumsum( )作出累积和; > y<-ifelse(x==0,-1,1 )
> cumsum(y)
(4)使用命令plot( ) 作出随机游动的示意图. > plot(cumsum(y))
3.在同一个图形中画出统计的四大分布密度曲线(dnorm, dchisq, dt, df)，注意不同分布有不
同的线型、颜色和宽度，还有图形都要在同一方框中，最后用图例说明(legend)。
> curve(dnorm(x,0,1),xlim=c(-1,5),ylim=c(0,0.5),col=1,lwd=1,lty=1)
> curve(dchisq(x,1),xlim=c(-1,5),ylim=c(0,0.5),lwd=2,lty=2,col=2,add=T)
> curve(dt(x,1),xlim=c(0,8),ylim=c(0,0.5),lwd=3,lty=3,col=3,add=T)
> curve(dt(x,1,1),xlim=c(0,8),ylim=c(0,0.5),lwd=4,lty=4,col=4,add=T)
> legend('topright',c("dnorm","dchisp","dt","df"),lty=c(1,2,3,4),col=c(1,2,3,4),lwd=c(1,2,3,4))
> curve(dnorm(x,0,1),xlim=c(-1,5),ylim=c(0,0.5),col=1,lwd=1,lty=1)
> curve(dchisq(x,1),xlim=c(-1,5),ylim=c(0,0.5),lwd=2,lty=8,col=2,add=T)
> curve(dt(x,1),xlim=c(0,8),ylim=c(0,0.5),lwd=5,lty=3,col=7,add=T)
> curve(dt(x,1,1),xlim=c(0,8),ylim=c(0,0.5),lwd=4,lty=4,col=4,add=T)
> legend('topright',c("dnorm","dchisp","dt","df"),lty=c(1,8,3,4),col=c(1,2,7,4),lwd=c(1,2,5,4))
>
4. 除本章给出的标准分布外, 非标准的随机变量X的抽样可通过格式点离散化方法实现.
设p (x )为X 的密度函数, 其抽样步骤如下
(1) 在X 的取值范围内等间隔地选取N 个点x 1, x 2,…, x N , 例如取N =1000; (2) 计算p (x i ); i = 1, 2, …, N ;
(3) 正则化p (x i ); i =1, 2,…,N , 使其成为离散的分布律, 即每一项除以
∑=N
i i
x p 1
)(;
(4) 按离散分布抽样方法使用命令sample( )从x i , i = 1, 2, … ,N 有放回地抽取n 个数, 例如 n =1000.
注：前面４小步是用来编一个函数，功能是对给定的概率密度产生随机数，形式应与ｒｎｏｒｍ差不多。
试以标准正态分布为例来说明. 为与R 中的正态抽样函数rnorm( )进行比较, 将作图区域分为左右两部分,
(1) 使用rnorm( )抽取n =1000个标准正态随机数, 并在左侧区域画出相应的直方图和核
密度估计曲线;
(2) 用格子点离散化抽样方法完成抽样, 并在右侧区域画出相应的直方图和核密度估计
曲线, 离散化所用的N =1000, n =1000, 取点范围为[-4, 4]. > N=seq(-4,4,length=1000) > f<-function(x)
+ dnorm(x)/sum(dnorm(x)) > f1=f(N)
> result=sample(N,replace=T,size=1000,prob=f1) > rn=rnorm(1000)
> op<-par(mfrow=c(1,2)) > hist(rn,probability=T)
> lines(density(rn),col="green",lwd=4) > hist(result,probability=T)
> lines(density(result),col="red",lwd=4) > par(op)。