R语言入门 S统计模拟.pptVIP

用S作随机模拟计算 作为统计工作者，我们除了可以用S迅速实 现新的统计方法，还可以用S进行随机模拟。 随机模拟可以验证我们的算法、比较不同 算法的优缺点、发现改进统计方法的方向， 是进行统计研究的最有力的计算工具之一。 下面先对模拟作一简单介绍 Buffon投针实验 模型：假设平面上有无数条距离为1的等距平行线,现向该平面随机地投掷长度ls=1的细针,则可以计算针与任一平行线相交的概率。 此处随机投针可以这样理解：针的中心与最近的平行线间的距离x均匀地分布在区间[0,1/2]上，针与平行线的夹角α均匀地分布在区间[0,π]上。于是针与线相交的充要条件是x/sin(α)=ls/2，从而针线相交的概率为 p=2*ls/π 随机模拟最基本的需要是产生伪随机数，S中已提供了大多数常用分布的伪随机数函数，可以返回一个伪随机数序列向量。 这些伪随机数函数以字母r开头，比如rnorm()是正态伪随机数函数，runif()是均匀分布伪随机数函数，其第一个自变量是伪随机数序列长度n。关于这些函数可以参见系统帮助文件。 下例产生1000个标准正态伪随机数： y - rnorm(1000) 这些伪随机数函数也可以指定与分布有关的参数，比如下例产生1000个均值为150、标准差为100的正态伪随机数： y - rnorm(1000, mean=150, sd=100) 产生伪随机数序列是不重复的，实际上，S在产生伪随机数时从一个种子出发，不断迭代更新种子，所以产生若干随机数后内部的随机数种子就已经改变了。有时我们需要模拟结果是可重复的，这只要我们保存当前的随机数种子，然后在每次产生伪随机数序列之前把随机数种子置为保存值即可： the.seed - .Random.seed …………… .Random.seed - the.seed y - rnorm(1000) 作为例子，我们来产生服从一个简单的线性回归 的数据。 # 简单线性回归的模拟 lm.simu - function(n){ # 先生成自变量。假设自变量x的取值范围在150到180之间，大致服从正态分布。 x - rnorm(n, mean=165, sd=7.5) # 再生成模型误差。假设误差服从N(0, 1.2) eps - rnorm(n, 0, 1.2) # 用模型生成因变量 y - 0.8 * x + eps return(data.frame(y,x)) } S没有提供多元随机变量的模拟程序，这里给出一个进行三元正态随机变量模拟的例子。假设要三元正态随机向量 的 n个独立观测，可以先产生n个服从三元标准正态分布的观测，放在一个 n行3列的矩阵中： U - matrix(rnorm(3*n), ncol=3, byrow=T) 可以认为矩阵U的每一行是一个标准的三元正态分布的观测。设矩阵 的Choleski分解 ， A为上三角矩阵，若随机向量 ，则 。因此， 作为一个三行 n列的矩阵每一行都是服从 分布的，且各行之间独立。经过转置，产生的 X X - matrix(rep(mu,n), ncol=3, byrow=T) + U %*% A 是一个 n行三列的矩阵。 有时模拟需要的计算量很大，多的时候甚至要计算几天的时间。对于这种问题我们要善于把问题拆分成可以单独计算的小问题，然后单独计算每个小问题，把结果保存在S对象中或文本文件中，最后综合保存的结果得到最终结果。 如果某一个问题需要的计算时间比较长，我们在编程时可以采用以下的技巧：每隔一定时间就显示一下任务的进度，以免计算已经出错或进入死循环还不知道；应该把中间结果每隔一段时间就记录到一个文本文件中（cat()函数可以带一个file参数和append参数，对这种记录方法提供了支持），如果需要中断程序，中间结果可能是有用的，有些情况下还可以根据记录的中间结果从程序中断的地方继续执行。 补充内容介绍 检验一个数据集的分布 两组数据同分布检验 检验一个数据集的分布 我们可以用很多方法分析一个单变量数据集的分布。最简单的办法就是直接看数字。利用函数summary 和fivenum 会得到两个稍稍有点差异的汇总信息。 此外，stem(茎叶图)也会反映整个数据集的数字信息。 attach(faithful) summary(eruptions) Min. 1st Qu. Median