统计模拟与R第一讲

比如计算积分
b
g( x)dx
a
其R程序如下：
也可采用R的内置函数来 计算积分
f1=function(n,a,b,g){
X=runif(n)
integrate(g,a,b)
sum((b-a)*g(a+(ba)*X))/length(X)
}
现举一个具体的例子比较他们的结果：
2
例 估计积分
e x x2 dx
plot(c(0,1,1,0),c(0,0,1,1),xlab=" ",ylab=" ") text(0,1,labels="A",adj=c(0,0)) text(1,1,labels="B",adj=c(1.5,0.5)) text(1,0,labels="C",adj=c(0.3,-0.8)) text(0,0,labels="D",adj=c(-0.5,0.1)) points(0.5,0.5) text(0.5,0.5,labels="o",adj=c(-1,0.3)) delta_t=0.01;n=200 x=matrix(0,nrow=5,ncol=n);x[,1]=c(0,1,1,0,0) y=matrix(0,nrow=5,ncol=n);y[,1]=c(1,1,0,0,1) d=c(0,0,0,0) for(j in 1:(n-1)){
解：建立平面直角坐标系，以时间间隔 t 进行采样，在每一时刻t计算每个人在下一
时刻 t t 时的坐标。不妨假设甲追踪乙，设在时刻t，甲、乙的坐标分别为
( x1 , y1 ), ( x2 , y2 ) 甲在 t t 时刻的坐标为 ( x1 vt cos , y1 vt sin ) d ( x2 x1 )2 ( y2 y1 )2 cos x2 x1 d , sin y2 y1 d ,
00
2.用模拟方法估计 Cov(U , eU ),其中U为(0,1)均匀分布的随机变量,并与其精确值
比较.
3.令 其中
xn 3xn1 5xn2 mod(100), n 3 x1 23, x2 66
求其前14个值.
例：追逐问题，在正方形ABCD的四个顶点各有一人， 在某一时刻，四人同时出发以匀速v走向顺时针方向的 下一个人。如果他们的方向始终保持对准目标，则最终 将按螺旋状曲线汇合于中心点O.试求这种情况下每个人的 轨迹。
}
例： 的估计:
也可以用平均值法:
1
0
1 x2 dx
4
程序
Pi2=function(n){ x=runif(n) 4*sum(sqrt(1-x^2))/n }
练习:1.估计下列积分
(1) 1 eex dx 0
(2) x(1 x2 )2dx 0
1
(3)
1 e( x y)2 dxdy
1, 0 x 1, 0 y 1 f ( x, y) 0, otherwise
则 P{ X 2 Y 2 1}
4
Pi<-function(n){ set.seed(1234579) k<-0; x<-runif(n); y<-runif(n); k <- length(x[x^2+y^2 < 1]) data.frame(Pi=4*k/n)
其中a和m为给定的正整数，每个 xn 均为0,1,2,…,m-1中的一个，
xn / m 称为一个伪随机数，它近似服从(0,1)上的均匀分布。
通常为了避免随机数的重复出现，可选取m为与计算机字长
相当的素数，如 m 231 1
பைடு நூலகம்
（2）混合同余法：利用如下递推公式
xn (axn1 c) mod m
例：对于(0,1)上均匀分布的随机变量 U1,U2,L , 定义
n
N min{n : Ui 1}
通过生成N来估计Ei[N1 ]
程序:f1.R
f1=function(n){ M=numeric(n)
练习:设 Ui , i 1 为随机数,定义N为
for(i in 1:n){ S=0;N=0 while(S<=1){
for(i in 1:4){ d[i]=sqrt((x[i+1,j]-x[i,j])^2+(y[i+1,j]-y[i,j])^2) x[i,j+1]=x[i,j]+delta_t*(x[i+1,j]-x[i,j])/d[i] y[i,j+1]=y[i,j]+delta_t*(y[i+1,j]-y[i,j])/d[i]
用混合同余法产生 xn 的R程序如下：
f=function(x0,a,c,m,n){ x=rep(0,n) x[1]=a*x0+c-m*floor((a*x0+c)/m for(i in 2:n){
x[i]=a*x[i-1]+c-m*floor((a*x[i-1]+c)/m) }
x }
二、利用随机数计算积分
第三章 随机数的产生及应用
一、乘同余法以及混合同余法产生随机数
随机数是指（0，1）上均匀分布随机变量的值，最早是手工 或通过机械方式由手纺车、投掷或洗纸牌等产生，现可由计 算机相继产生伪随机数 （1）乘同余法：一种通常的产生伪随机数的方法为：
从初值或者称为种子 x0 出发，利用公式
xn axn1 mod m
n
N max{n : Ui e3 } i 1
S=S+runif(1) N=N+1
用模拟法求E[N]
}
M[i]=N }
用模拟法求P{N i}, i 0,1, 2, 3,4,5,6
mean(M)
}
例： 的估计:
考虑服从(0,1)区间上均匀分布的独立的随机变量
X ,Y 因此,二维随机变量(X,Y)的联合概率密度为
随机数的最早应用之一就是积分计算。设g(x)是一个函数，计算
1
g(x)dx
0
注意到若U在(0,1)上服从均匀分布，则
E[g(U )]
产生k个随机数 U1,U2 ,L ,Uk 利用强大数定律，我们以概率1有
k g(Ui ) E[g(U )] , k
i1 k
近似计算上述积分的方法称为Monte Carlo方法。
2
> f1(100000,-2,2,g) [1] 92.05593
> integrate(g,-2,2) 93.16275 with absolute error < 0.00062
11
例：估计二重积分 e(x y)2dxdy
00
给出如下算法： X=runif(100000,0,1) Y= runif(100000,0,1) f=function(x,y){exp((x+y)^2)} sum(f(X,Y))/length(X) [1] 4.900897