R语言入门 S统计模拟

1、实验目的实验目的 实验内容实验内容 学习如何应用学习如何应用R软件进行简单统计模拟软件进行简单统计模拟 1 1、统计模拟简介、统计模拟简介 2 2、应用实例、应用实例 S统计分析统计分析 用用S作随机模拟计算作随机模拟计算 作为统计工作者，我们除了可以用S迅速实 现新的统计方法，还可以用S进行随机模拟。 随机模拟可以验证我们的算法、比较不同 算法的优缺点、发现改进统计方法的方向， 是进行统计研究的最有力的计算工具之一。 下面先对模拟作一简单介绍 模拟的概念模拟的概念 模拟就是利用物理的、数学的模型来类比、模 仿现实系统及其演变过程，以寻求过程规律的一 种方法。 模拟的基本思想是建立一个试验模型

2、，这个模 型包含所研究系统的主要特点通过对这个实验 模型的运行，获得所要研究系统的必要信息。 模拟的方法模拟的方法 1、物理模拟物理模拟： 对实际系统及其过程用功能相似的实物系统 去模仿。例如，军事演习、船艇实验、沙盘 作业等。 物理模拟通常花费较大、周期较长，且在 物理模型上改变系统结构和系数都较困难。 而且，许多系统无法进行物理模拟，如社会 经济系统、生态系统等。 在实际问题中，面对一些带随机因素的复杂系 统，用分析方法建模常常需要作许多简化假设， 与面临的实际问题可能相差甚远，以致解答根本 无法应用。这时，计算机模拟几乎成为唯一的选 择。 在一定的假设条件下，运用数学运算模拟系统在一定的

3、假设条件下，运用数学运算模拟系统 的运行，称为数学模拟。现代的数学模拟都是在的运行，称为数学模拟。现代的数学模拟都是在 计算机上进行的，称为计算机模拟。计算机上进行的，称为计算机模拟。 2、数学模拟数学模拟 计算机模拟可以反复进行，改变系统的结构和系 数都比较容易。 随机模拟方法随机模拟方法是一种应用随机数来 进行计算机模拟的方法此方法对研究 的系统进行随机抽样观察，通过对样本 值的观察统计，求得所研究系统的某些 参数 随机模拟方法随机模拟方法也称为蒙特卡洛方法,其 起源最早可以追溯到18世纪下半叶的Buffon 试验. Buffon投针实验 模型：模型：假设平面上有无数条距离为1的等 距平行

4、线,现向该平面随机地投掷长度 ls=1的细针,则可以计算针与任一平行线 相交的概率。 此处随机投针可以这样理解：针的中心与 最近的平行线间的距离x均匀地分布在区 间0,1/2上，针与平行线的夹角均匀地分 布在区间0,上。于是针与线相交的充要 条件是x/sin() y y the.seed .Random.seed y - rnorm(1000) 作为例子，我们来产生服从一个简单的线性回归 的数据。 # 简单线性回归的模拟 lm.simu - function(n) # 先生成自变量。假设自变量x的取值范围在 150到180之间，大致服从正态分布。 x - rnorm(n, mean=165,

5、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为上三角

6、矩阵，若随机向量 ， 则 。因此， 作为一个三行 n列的矩阵每一行都是服从 分布的，且各行之间独立。经过转置，产生的 X ( , )XN T A A 3 (0,)NI ( , ) T AN ( , )N ( , ) T A UN X attach(faithful) summary(eruptions) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.600 2.163 4.000 3.488 4.454 5.100 fivenum(eruptions) 1 1.6000 2.1585 4.0000 4.4585 5.1000 stem(eruptions) 茎

7、叶图和柱状图相似，R 用函数hist 绘制柱状图。 hist(eruptions) # 让箱距缩小，绘制密度图 hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE) lines(density(eruptions, bw=0.1) rug(eruptions) # 显示实际的数据点 更为精致的密度图是用函数density 绘制的。在这个例子 中，我们加了一条由density 产生的曲线。 你可以用试错法（trial-and-error）选择带宽bw （bandwidth）。因为默认的带宽值让密度曲线过于平滑。 (现在已经有一些自动的带宽挑选方法) 我们可

8、以用函数ecdf 绘制一个数据集的经验累积分布 （empirical cumulative distribution）函数。 plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE) 显然，这个分布和其他标准分布差异很大。那么右边的 情况怎么样呢，就是火山爆发3分钟后的状况？我们可 以拟合一个正态分布，并且重叠前面得到的经验累积密 度分布。 long 3 plot(ecdf(long), do.points=FALSE, verticals=TRUE) x lines(x, pnorm(x, mean=mean(long), sd=sqrt(

9、var(long), lty=3) 分位比较图（Quantile-quantile (Q-Q) plot）便于我们更 细致地研究二者的吻合程度。 par(pty=s) # 设置一个方形的图形区域 qqnorm(long); qqline(long) 上述命令得到的QQ图表明二者还是比较吻合的， 但右侧尾部偏离期望的正态分布 我们可以用t 分布获得一些模拟数据以重复上面的 过程 x - rt(250, df = 5) qqnorm(x); qqline(x) 这里得到的QQ图常常会出现偏离正态期望的长尾 区域(如果是随机样本)。 我们可以用下面的命令针对特定的分布绘制Q-Q图 qqplot(qt

10、(ppoints(250), df = 5), x, xlab = Q-Q plot for t dsn) qqline(x) 最后，我们可能需要一个比较正规的正态性检 验方法。 R提供了Shapiro-Wilk 检验 shapiro.test(long) 和Kolmogorov-Smirnov 检验 ks.test(long, pnorm, mean = mean(long), sd = sqrt(var(long) 注意一般的统计分布理论（distribution theory） 在这里可能无效，因为我们用同样的样本对正 态分布的参数进行估计的。 双样本检验双样本检验 前面已经学会了单样本

11、的正态性检验。而更 常见的是比较两个样本的特征。在R 里面， 所有“传统”的检验都放在包stats 里面。这 个包常常会自动载入。 下面是冰融化过程的潜热（latent heat） (cal/gm) 数据（来自Rice (1995, p.490)） Method A: 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97 80.05 80.03 80.02 80.00 80.02 Method B: 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97 盒状图（boxplot）为这两组数据提供了简单 的图形比较

12、。 A - scan() 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97 80.05 80.03 80.02 80.00 80.02 B t.test(A, B) 上面的结果表明在正态前提下，二者有明 显的统计差异。R 函数默认两个样本方差 不齐，而SPLUS 类似函数t.test 则默认方差 齐性。如果两个样本都是来自正态总体， 我们可以用F检验来确定方差的齐性情况， var.test(A, B) 这表明二者方差在统计学上没有显著差异， 我们可以采用传统的假设方差齐性的t-检验。 注:t-检验有正态性假设的, 所以在进行t-检 验前, 原则上先要做数据的正态性检验. 所有这些检验都假设了数据的正态性。 双样本的Wilcoxon (或者Mann-Whitney) 检验没有正态性的前提，仅仅要求在原假 设(null hypothesis)情况下样本来自一个常 规的连续分布。 wilcox.test(A, B) 注意警告信息：在两个样本中都有同秩现 象, 这表明这些数据来自离散分布(可能由于 数据的近似处理造成)。 有好多种方法可以图形化的显示两个样本 的差别。我们已经看过盒状图的比较。 下面的命令 plot(ecdf(A), do.points=FALSE, verticals=TRUE, xlim=range(A, B)