R语言实验报告6 - 图文

一、 实验目的

1. 用 R 生成服从某些具体已知分布的随机变量

二、 实验内容

在 R 中各种概率函数都有统一的形式，即一套统一的前缀+分布函名：

d 表示密度函数（density）；

p 表示分布函数（生成相应分布的累积概率密度函数）；

q 表示分位数函数，能够返回特定分布的分位（quantile）； r 表示随机函数，生成特定分布的随机数（random）。

1、 通过均匀分布随机数生成概率分布随机数的方法称为逆变换法。对于任意随机变量 X，其分布函数为 F，定义其广义逆为：F-(u)=inf{x;F(x)≥u}若u~u (0,1)，则F-（u）和X 的分布一样

Example 1 如果X~Exp（1）（服从参数为 1 的指数分布），F（x）=1-e-x。若u=1-e-x并且u~u(0,1)，则X=-logU~Exp(1)

则可以解出x=-log(1-u)

通过随机数生成产生的分布与本身的指数分布结果相一致 R 代码如下：

nsim = 10^4

X = -log(U) U = runif(nsim)

Y = rexp(nsim) X11(h=3.5)X

par(mfrow=c(1,2),mar=c(2,2,2,2))

curve(dexp(x),add=T,col=\

hist(X,freq=F,main=\

hist(Y,freq=F,main=\curve(dexp(x),add=T,col=\

2、 某些随机变量可由指数分布生成。若X i~ Exp(1) 独立同分布的随机变量，那么从X i出

发可以得到以下三个标准分布

Example 2 生成自由度为 6 的χ2分布。

R 代码如下：

nsim = 10^4

U=matrix(data=U,nrow=3) X=-log(U)

X=2* apply(X,2,sum) Y = rchisq(nsim,df = 6)

3、 一种正态分布随机变量模拟使用 Box-Muller 算法得到 N(0,1)随机变量。这种方法与基于

中心极限定理的近似算法相比较，Box-Muller 算法是精确的，它由两个均匀分布产生两个独立的正态分布，其仅有的缺点是必须计算 log，cos，sin。具体解释如下：如果两个随机变量U1, U2独立同分布于u(0,1)，那么由此可以产生两个独立的正态分布X1, X2。

par(mfrow=c(1,2),mar=c(2,2,2,2)) hist(X,freq=F,main=\

Exp\plot(d)

Example 3 两次产生相同的 10000 个服从正态分布的随机数，作其中一个的概率直方图，并添加正态分布的密度函数线。

nsim = 10^4

set.seed(22)

x1 <- rnorm(nsim,mean = 0, sd = 1) set.seed(22)

x2 <- rnorm(nsim,mean = 0, sd = 1)

这里的 x1，x2 数值完全相同。

4、 设随机变量 X 的分布列P{X=xi}=pi, 记p(0)=P{X≤0}=0, … , p(i)=P{X≤xi}=

?ij?1pi。

设 r 是[0,1]区间上的均匀分布的随机数。当且仅当p(i-1)

Example 4 生成二项分布的随机数

首先生成二点分布 > size = 1;p = 0.5

> rbinom(10,size,p) [1] 1 0 0 0 0 1 1 0 1 0

接下来生成服从 B（10，0.5）的二项分布 > size = 10; > p = 0.5

> rbinom(5,size,p) [1] 6 6 5 8 5

由此可见，随着实验次数 n 的增大，二项分布越来越接近正态分布。 R代码如下 size = 1;p = 0.5 rbinom(10,size,p) #生成 5 个服从 B（10，0.5）的二项分布随机数 size = 10; p = 0.5 rbinom(5,size,p) par(mfrow = c(1,3)) p=0.25 for(n in c(10,20,50)) {x <- rbinom(100,n,p) hist(x,prob = T,main = paste(\ xvals = 0:n points(xvals,dbinom(xvals,n,p),type = \ } par(mfrow = c(1,1))

5．服从 beta 分布的随机变量的一般算法可以基于 Accept-Reject method 来生成，

用的 工具分布为均匀分布 U[0,1]，假设两个参数都大于 1（通用的 rbeta 函数无此约束）。上 界 M 为 beta 分布密度的最大值， Example 5 对于α =2.7，β =6.3 ，最大值 M=2.67，求出 beta 分布 R代码 Nsim=2500 a=2.7;b=6.3 M=2.67 u=runif(Nsim,max=M) y=runif(Nsim) x=y[u

Accept-Reject algorithm的几个关键要点：

1. 仅要求比率f/M，所有算法并不依赖于正则化常数。

2. f/Mg的上界不必是紧的，该算法一人有效，即M可由任何更大的常数替代。

3.接受概率为1/M，故M应该尽可能小。

Accept-Reject algorithm 的效率可由接受概率度量，此概率越高，由g模拟的数据浪费的越少。这当然也必须考虑从g生成模拟值的计算成本。故要选择对f最好的g。