摘要

除了概率逆变换，还有多种变换可以用来生成随机变量的样本。接下来会给出一些变换以及两个示例。

文章目录

  • 摘要
  • 几个变换
  • 实例
    • Beta分布
    • 例2 对数分布

几个变换

  1. 如果 Z ∼ N ( 0 , 1 ) Z\sim N(0,1) Z∼N(0,1)，则 V = Z 2 ∼ χ 2 ( 1 ) V = Z^2 \sim\chi ^2(1) V=Z2∼χ2(1)
  2. 如果 U ∼ χ 2 ( m ) U \sim \chi^2(m) U∼χ2(m)与 V ∼ χ 2 ( n ) V \sim \chi^2(n) V∼χ2(n)相互独立，则 F = U / m V / n F = \frac{U/m}{V/n} F=V/nU/m​服从自由度为(m,n)的F分布
  3. 如果 Z ∼ N ( 0 , 1 ) Z\sim N(0,1) Z∼N(0,1)与 V ∼ χ 2 ( n ) V \sim \chi^2(n) V∼χ2(n)相互独立，则 T = z V / n T= \frac{z}{\sqrt{V/n}} T=V/n ​z​服从自由度为 n 的F分布
  4. 如果 U , V ∼ U n i f ( 0 , 1 ) U,V\sim Unif(0,1) U,V∼Unif(0,1)且相互独立，则
    Z 1 = − 2 log ⁡ U cos ⁡ 2 π V Z 2 = − 2 log ⁡ V cos ⁡ 2 π U Z_1 = \sqrt{-2\log{U}}\cos{2\pi V} \\ Z_2 = \sqrt{-2\log{V}}\cos{2\pi U} Z1​=−2logU ​cos2πVZ2​=−2logV ​cos2πU
    是独立的标准正态分布
  5. 如果 U ∼ G a m m a ( r , λ ) 与 V ∼ G a m m a ( s , λ ) U\sim Gamma(r,\lambda)与V\sim Gamma(s,\lambda) U∼Gamma(r,λ)与V∼Gamma(s,λ)相互独立，则 X = U U + V X = \frac{U}{U+V} X=U+VU​服从beta(r, s)分布。
  6. 如果 U , V ∼ U n i f ( 0 , 1 ) U,V\sim Unif(0,1) U,V∼Unif(0,1)且相互独立，则
    X = [ 1 + log ⁡ V log ⁡ ( 1 − ( 1 − θ ) 2 ) U ] X=[1+\frac{\log{V}}{\log({1-(1-\theta)^2)^U}}] X=[1+log(1−(1−θ)2)UlogV​]
    服从对数分布 log ⁡ ( θ ) , [ X ] \log (\theta),[X] log(θ),[X]表示X的整数部分。

基于变化5，6的生成器会举例演示。另外，求和与混合作为两种特殊的变换，会另外单独说明。

实例

Beta分布

下面有关Beta与Gamma分布的关系式提供了另一种方式的随机数生成器：

如果 U ∼ G a m m a ( r , λ ) 与 V ∼ G a m m a ( s , λ ) U\sim Gamma(r,\lambda)与V\sim Gamma(s,\lambda) U∼Gamma(r,λ)与V∼Gamma(s,λ)相互独立，则 X = U U + V X = \frac{U}{U+V} X=U+VU​服从beta(r, s)分布。

这个转换决定了一个生成服从Beta(a,b)的随机样本算法，具体步骤如下：

  1. 生成一个服从Gamma(a,1)的随机样本
  2. 生成来自分布gamma(b,1)的随机样本v
  3. 传递赋值 x = u u + v x = \frac{u}{u+v} x=u+vu​

代码实现：

n <- 1000 
a <- 3 
b <- 2 
u <- rgamma(n, shape=a, rate=1) 
v <- rgamma(n, shape=b, rate=1) 
x <- u / (u + v)

该样本数据可以利用一个QQ分位数图与Beta(3,2)分布进行比较。如果采样分布是Beta(3,2)，QQ图应该是接近线性的。

q <- qbeta(ppoints(n), a, b) 
qqplot(q, x, cex=0.25, xlab="Beta(3, 2)", ylab="Sample") 
abline(0, 1)

添加线x=q供参考。

图3.2中有序样本与Beta(3,2)分位数QQ图是非常接近线性的，因为如果生成的样本实际上是一个Beta(3,2)样本，那么它应该和它差不多。

例2 对数分布

这个例子提供了另一个更有效率的对数分布随机数的生成器。(前面逆变换法中也有一个生成服从对数分布随机样本的方法，比这个麻烦一点，后面再补上。)
如果 U , V ∼ U n i f ( 0 , 1 ) U,V\sim Unif(0,1) U,V∼Unif(0,1) 且相互独立，则
X = [ 1 + log ⁡ V log ⁡ ( 1 − ( 1 − θ ) 2 ) U ] X=[1+\frac{\log{V}}{\log({1-(1-\theta)^2)^U}}] X=[1+log(1−(1−θ)2)UlogV​]
服从对数分布log(θ), 这个变换提供了一个简单有效的对数分布随机数生成办法：

  1. 生成服从U(0,1)分布的样本u
  2. 生成服从U(0,1)分布的样本v
  3. 传递赋值 x = [ 1 + log ⁡ v log ⁡ ( 1 − ( 1 − θ ) 2 ) u ] x=[1+\frac{\log{v}}{\log({1-(1-\theta)^2)^u}}] x=[1+log(1−(1−θ)2)ulogv​]

以下是对数分布log(0.5)与利用转换方法生成的样本的比较。

− 1 l o g ( 1 − θ ) ∗ θ k / k \frac{-1}{log(1 - \theta)} * \theta^k / k log(1−θ)−1​∗θk/k

n <- 1000 
theta <- 0.5 
u <- runif(n) #第1步
v <- runif(n) #第2步
x <- floor(1 + log(v) / log(1 - (1 - theta)^u)) #第3步
k <- 1:max(x)  #calc. logarithmic probs. 
p <- -1 / log(1 - theta) * theta^k / k  
se <- sqrt(p*(1-p)/n) 
p.hat <- tabulate(x)/n  
#经验概率p.hat在理论概率p的两个标准误差范围内。
> print(round(rbind(p.hat, p, se), 3))
       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]
p.hat 0.728 0.183 0.058 0.021 0.005 0.004 0.001
p     0.721 0.180 0.060 0.023 0.009 0.004 0.002
se    0.014 0.012 0.008 0.005 0.003 0.002 0.001

tabulate函数为正整数，所以它可以在对数样本上使用。请将数据重新编码为正整数或使用table。如果数据不是正整数，tabulate将截断实数，并无警告地忽略小于1的整数。

下面这个函数是一个简单的对数函数

rlogarithmic <- function(n, theta) { 
	stopifnot(all(theta > 0 & theta < 1))   
	th <- rep(theta, length=n) 
	u <- runif(n) 
	v <- runif(n) 
	x <- floor(1 + log(v) / log(1 - (1 - th)^u)) 
	return(x) 
}

更多推荐

分布式统计计算-----变换法生成随机数R studio