如何使用Box-Muller方法生成卡方与t分布随机数及R代码实现

卡方分布生成问题解答 #

你可以直接将Box-Muller生成的X1、X2分别平方后相加得到卡方分布随机值，原理如下：
Box-Muller方法输出的X1、X2是互相独立的*标准正态分布（N(0,1)）*随机变量，而自由度为k的卡方分布的定义就是k个独立标准正态变量的平方和服从的分布，因此X1^2 + X2^2服从自由度为2的卡方分布。
如果需要生成更高自由度的卡方分布，只需要生成对应数量的独立标准正态变量，将所有变量平方求和即可。

示例代码：生成1000个自由度为2的卡方随机值

n <- 1000
A <- runif(n, 0, 1)
B <- runif(n, 0, 1)
X1 <- sin(2*pi*A)*sqrt(-2*log(B))
X2 <- cos(2*pi*A)*sqrt(-2*log(B))
# 自由度为2的卡方分布随机值
chi2_df2 <- X1^2 + X2^2

t分布生成方法 #

自由度为k的t分布定义为：若Z ~ N(0,1)（标准正态分布）、V ~ χ²(k)（自由度为k的卡方分布），且Z和V互相独立，则变量T = Z / sqrt(V/k)服从自由度为k的t分布。
基于已有的Box-Muller实现，生成逻辑如下：

• 生成独立的标准正态变量作为分子Z
• 生成k个独立标准正态变量，平方求和得到自由度为k的卡方变量V，注意V需要和Z互相独立
• 按公式计算得到t分布随机值

示例代码：生成1000个自由度为3的t分布随机值

n <- 1000
# 生成标准正态变量Z（分子）
A_z <- runif(n, 0, 1)
B_z <- runif(n, 0, 1)
Z <- sin(2*pi*A_z)*sqrt(-2*log(B_z))
# 生成3个独立标准正态，计算自由度为3的卡方变量V（和Z独立）
A1 <- runif(n, 0, 1)
B1 <- runif(n, 0, 1)
norm1 <- sin(2*pi*A1)*sqrt(-2*log(B1))
norm2 <- cos(2*pi*A1)*sqrt(-2*log(B1))
A2 <- runif(n, 0, 1)
B2 <- runif(n, 0, 1)
norm3 <- sin(2*pi*A2)*sqrt(-2*log(B2))
V <- norm1^2 + norm2^2 + norm3^2
# 自由度为3的t分布随机值
t_df3 <- Z / sqrt(V/3)

注意事项 #

计算卡方变量V使用的正态随机值，不能和分子Z使用的正态随机值来自同一批Box-Muller输出，否则会违反独立性假设，导致生成的随机值不符合预期分布。

内容的提问来源于stack exchange，提问作者George