R语言实验报告

一、实验目的
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
U = runif(nsim)
X = -log(U)
Y = rexp(nsim)
X11(h=3.5)X
par(mfrow=c(1,2),mar=c(2,2,2,2))
hist(X,freq=F,main="Exp from Uniform",ylab="",xlab="",ncl=150,col="grey",xlim=c(0,8))
curve(dexp(x),add=T,col="sienna",lwd=2)
hist(Y,freq=F,main="Exp from R",ylab="",xlab="",ncl=150,col="grey",xlim=c(0,8))
curve(dexp(x),add=T,col="sienna",lwd=2)
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)
par(mfrow=c(1,2),mar=c(2,2,2,2))
hist(X,freq=F,main="chisq from
Exp",ylab="",xlab="",ncl=150,col="grey",xlim=c(0,8)) d = density(Y)
plot(d)
3、一种正态分布随机变量模拟使用Box-Muller 算法得到N(0,1)随机变量。

这种方法与基于
中心极限定理的近似算法相比较，Box-Muller算法是精确的，它由两个均匀分布产生两个独立的正态分布，其仅有的缺点是必须计算log，cos，sin。

具体解释如下：如果两个随机变量U1, U2独立同分布于u(0,1)，那么由此可以产生两个独立的正态分布X1, X2。

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=x i}=p i, 记p(0)=P{X≤0}=0, … , p(i)=P{X≤x i}=∑=i j i p1。

设 r 是[0,1]区间上的均匀分布的随机数。

当且仅当p(i-1)<r<p(i)时，令X=x i，则P{p(i-1)<r<p(i)}=p(i)-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("n = ",n))
xvals = 0:n
points(xvals,dbinom(xvals,n,p),type = "h",lwd = 3)
}
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<dbeta(y,a,b)]
这里的 x 服从 beta 分布
Accept-Reject algorithm的几个关键要点：
1. 仅要求比率f/M，所有算法并不依赖于正则化常数。

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

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

Accept-Reject algorithm 的效率可由接受概率度量，此概率越高，由g模拟的数据浪费的越少。

这当然也必须考虑从g生成模拟值的计算成本。

故要选择对f最好的g。