R软件 蒙特卡罗模拟.doc

R使用指南打开R下图是R软件的主窗口，R软件的界面与Windows的其他编程软件类似，由一些菜单和快捷按钮组成。快捷按钮下面的窗口便是命令输入窗口,它也是部分运算结果的输出窗口,有些运算结果则会在新建的窗口中输出。当一个R 程序需要你输入命令时,它会显示命令提示符。默认的提示符是。技术上来说， R 是一种语法非常简单的表达式语言（expression language）。它大小写敏感，因此A 和a 是不同的符号且指向不同的变量。可以在R 环境下使用的命名字符集依赖于R 所运行的系统和国家(就是系统的locale 设置)。通常，数字，字母，. 和都是允许的(在一些国家还包括重音字母)。不过，一个命名必须以. 或者字母开头，并且以. 开头时第二个字符不允许是数字。基本命令要么是表达式（expressions）要么就是赋值（assignments）。如果一条命令是表达式，那么它将会被解析（evaluate），并将结果显示在屏幕上，同时清空该命令所占内存。赋值同样会解析表达式并且把值传给变量但结果不会自动显示在屏幕上。命令可以被(;)隔开，或者另起一行。基本命令可以通过大括弧(f和g) 放在一起构成一个复合表达式（compound expression）。注释几乎可以放在任何地方7。一行中，从井号(#)开始到句子收尾之间的语句就是注释。如果一条命令在一行结束的时候在语法上还不完整， R 会给出一个不同的提示符，默认是+。该提示符会出现在第二行和随后的行中，它持续等待输入直到一条命令在语法上是完整的。该提示符可以被用户修改。在后面的文档中，我们常常省略延续提示符（continuation prompt），以简单的缩进表示这种延续。R的帮助首先先来看看如何使用帮助文件 这里有两个方式：在R中输入help.start()通过启动HTML形式的在线帮助(使用你的计算机里面可用的浏览器)。你可以用鼠标点击上面的链接。如图：英文足够好的话，这就是R的最佳教材（可惜我不行！）。当你看到某个不知道语句可以通过这方式来找到答案：例如plot（）语句，你想知道他是干什么的那么就在R中输入“?plot()”看到了吗，很详细的解释。看懂了plot（）是做什么的了吗？，没看懂看下面的例子。下面的这部分会话让你在操作中对R 环境的一些特性有个简单的了解。你对系统的许多特性开始时可能有点不熟悉和困惑，但这些迷惑会很快消失的。登录，启动你的桌面系统。R 程序开始，并且有一段引导语。(在R 里面，左边的提示符将不会被显示防止混淆。)help.start()启动HTML 形式的在线帮助(使用你的计算机里面可用的浏览器)。你可以用鼠标点击上面的链接。最小化帮助窗口，进入下一部分。x - rnorm(50)y - rnorm(x)这里x-rnorm(50)产生50 个标准正态数据，y-rnorm(x)是产生和x为数一样多的一组标准正态数据。plot(x, y)画二维散点图。一个图形窗口会自动出现。ls()查看当前工作空间里面的R 对象。rm(x, y)去掉不再需要的对象。(清空)。x - 1:20等价于x = (1; 2; : : : ; 20)。w - 1 + sqrt(x)/2标准差的权重向量。dummy - data.frame(x=x, y= x + rnorm(x)*w)dummy创建一个由x 和y构成的双列数据框，查看它们。fm - lm(y x, data=dummy)summary(fm)拟合y 对x 的简单线性回归，查看分析结果。fm1 - lm(y x, data=dummy, weight=1/w2)summary(fm1)现在我们已经知道标准差，做一个加权回归。attach(dummy)让数据框中的列项可以像一般的变量那样使用。lrf name - function(arg 1 , arg 2 , .) expression其中expression 是一个R 表达式(常常是一个成组表达式)，它利用参数argi 计算最终的结果。该表达式的值就是函数的返回值。可以在任何地方以name(expr1 , expr2 , .) 的形式调用函数。计算机模拟方法Monte Carlo方法 下面主要介绍最基本的计算机模拟方法Monte Carlo方法,此方法的基本思想是将各种随机时间的概率特征(概率分布、数学期望)与随机事件的模拟联系起来,用试验的方法确定事件的相应概率与数学期望.因而, Monte Carlo方法的突出特点是概率模型的解是由试验得到的,而不是计算出来的。此外,模拟任何一个实际过程, Monte Carlo方法都需要用到大量的随机数,计算量很大,人工计算是不可能的,只能在计算机上实现.在概率论中,著名的Buffon掷针问题就是用统计试验的方法求圆周率的典型代表.现在用Monte Carlo方法模拟一下Buffon投针试验，对于该问题大家都已经了解了其解法,下面就不再赘述,我们只从随机模拟的角度度该问题进行求解。针与平行线相交的充分必要条件是。Buffon的投针试验在计算机上实现,需要一下两个步骤：产生随机数.首先产生n个相互独立的随机变量,x的抽样序列,其中.模拟检验：检验不等式是否成立.若上式成立,表示第i次试验成功(即针与平行线相交).设n次试验中有k次成功,则的估值为，其中,均为预先给定。将上述步骤编写成R模拟程序。模拟代码如下：buffon-function(n, l=0.8, a=1) #定义一个新的函数buffon,包括三个自变量(n ,l, a) k-0 #给k赋初值为0 theta-runif(n,0, pi); x-runif(n,0, a/2) #定义变量theta和x分别为(0,pi)和(0,1)之间的均匀分布的随机数 for (i in 1:n) #控制循环语句i从1到n if (xi= l/2*sin(thetai) #条件语句如果括号内的条件成立,继续做下一步 k buffon(100) 1 3.137255 buffon(1000) 1 3.137255 buffon(10000) 1 3.142801正态分布随机数的产生；这里主要介绍极限近似法:设是(0,1)区间上n个独立的均匀分布的随机数,由中心极限定理得到,近似地服从正态分布N(0,1),为了保证一定的精度,上式中的n应取得足够大,一般大约取n=10左右,为方便起见,可取n=12.此时,上式由最简单的形式。当是(0,1)上的随机数时,则也是(0,1)上的随机数,因此上式可以改写为。若随机数x服从N(0,1)分布时,令,则y是正态分布的随机数.由此可以得到任意参数的正态分布的随机数。用Monte Carlo方法模拟上述过程：ztsj=function(n, miu,sigma) #定义一个新的函数ztsj,包括三个自变量(n , miu,sigma) k=1;y=rep(0,n) #给k赋初值,定义y为一个n维0向量, for(k in 1:n) #控制循环语句i从1到n x=c(runif(12) #定义x为一个12维的行向量,每个元素都是均匀分布随机数yk=sigma*(sum(x)-6)+miu #计算向量y的每个分量k=k+1 #k进行运算k=k+1 y #输出y运行上述程序之后就在R的主窗口输入:ztsj(100,1,2)就会得到如下图的输出结果。当然,括号中的数字可以任意改写。用R软件生成随机数的方法：（1） runif产生均匀分布的随机数，参数为n，a，b,其中n为随机数的个数，a，b为区间（a，b）端点值，当a,b默认时,为(0,1)区间上的随机数。（2） rnorm产生正态分布的随机数,参数为n,其中n为随机数的个数, 为均值,为标准差,当,默认时,为标准正态分布N(0,1)的随机数。（3） rpois产生Poisson分布的随机数,参数为n,其中n为随机数的个数, 为Poisson分布的参数。R软件还可以产生其他分布的随机数,这里就不一一列举了,这些分布的函数前加r,就表示是生成该分布的随机数。在R中使用包（packages）在本节，我们将以mcmc过程为例，演示如何在R中使用第三方软件包。这里使用mcmc package ，将该压缩包安装到r中并加载程序包（程序包-加载程序包）用以