题目如下，总的来说就是一个标准正态随机向量的生成，并限制在近似抽样法上：
题目分析：
由于目标仅是为了产生50个独立同\(N(0,1)\)分布的随机变量，所以并不需要使用生成随机向量的方法来考虑变量之间的非独立性，只需要通过指定的方法生成N(0,1)的单个变量并循环50次即可。为了更好地展现每个方法的运行实现和生成随机数的效果，我们不妨合理修改，将题目要求的生成“50”个增加到“500”个。在指定的“近似抽样”方法中，我们可以使用一系列近似的构造方法，包括：
比较思路框架：
step1：在调用这些方法之前，我们需要为每个方法设定初始化参数，如题目中提示，然将这些初始化参数对应的基础分布和生成的正态分布随机数图像画出，用不同的检验方法来验证随机数生成的效果。
step2：灵敏度分析，这些方法相对他们各自的参数的灵敏度如何？我们首先确立指标，然后将更改不同的参数并记录相应的指标变化，最后进行比较。
step3：对于不同的方法，哪些算法在生成题目要求的随机向量时更有效率？我们将从时间和空间复杂度上进行比较。
基础分布的频数图绘制函数
library(ggplot2)
plot_basefun <- function(samples,name="",color="#ffc556") {
ggplot(data.frame(samples), aes(x = samples)) +
geom_histogram(aes(y = after_stat(density)), bins = 30, fill = color, alpha = 0.7) +
geom_density(color = "#abddff", linewidth = 1) +
labs(x = paste0(name,"随机数"), y = "概率密度") +
theme_minimal()}get_norm_by_exp <- function(lambda = 1, n = 12, m = 500) {
exp_samples <- replicate(m, sum(rexp(n, rate = lambda)))
norm_samples <- (exp_samples - n / lambda) / sqrt(n / lambda^2)
return(norm_samples)
}
plot_basefun(rexp(500, rate = 1),name="指数分布")
get_norm_by_unif <- function(min = 0, max = 1, n = 12, m = 500) {
unif_samples <- replicate(m, sum(runif(n, min = min, max = max)))
norm_samples <- (unif_samples - n * (min + max) / 2) / sqrt(n * (max - min)^2 / 12)
return(norm_samples)
}
plot_basefun(runif(5000, min = 0, max = 1),name="平均分布")
get_norm_by_gamma <- function(shape = 2, scale = 1, n = 12, m = 500) {
gamma_samples <- replicate(m, sum(rgamma(n, shape = shape, scale = scale)))
norm_samples <- (gamma_samples - n * shape * scale) / sqrt(n * shape * scale^2)
return(norm_samples)
}
plot_basefun(rgamma(500, shape = 2, scale = 1),name="伽马分布")
get_norm_by_beta <- function(shape1 = 2, shape2 = 5, n = 12, m = 500) {
beta_samples <- replicate(m, sum(rbeta(n, shape1 = shape1, shape2 = shape2)))
mean_beta <- shape1 / (shape1 + shape2)
var_beta <- (shape1 * shape2) / ((shape1 + shape2)^2 * (shape1 + shape2 + 1))
norm_samples <- (beta_samples - n * mean_beta) / sqrt(n * var_beta)
return(norm_samples)
}
plot_basefun(rbeta(500, shape1 = 2, shape2 = 5),name="贝塔分布")
library(ggplot2)
library(qqplotr)
Attaching package: 'qqplotr'The following objects are masked from 'package:ggplot2':
stat_qq_line, StatQqLinelibrary(gridExtra)
library(grid)
plot_distribution <- function(samples,name="") {
# 绘制频率直方图和密度曲线
hist_plot <- ggplot(data.frame(samples), aes(x = samples)) +
geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "#ffc556", alpha = 0.7) +
geom_density(color = "#abddff", linewidth = 1) +
labs(x = "随机数", y = "概率密度") +
theme_minimal()
# 绘制QQ图
qq_plot <- ggplot(data = data.frame(samples), mapping = aes(sample = samples)) +
geom_qq_band(bandType = "ks", mapping = aes(fill = "KS"), alpha = 0.5) +
geom_qq_band(bandType = "ts", mapping = aes(fill = "TS"), alpha = 0.5) +
geom_qq_band(bandType = "pointwise", mapping = aes(fill = "Normal"), alpha = 0.5) +
geom_qq_band(bandType = "boot", mapping = aes(fill = "Bootstrap"), alpha = 0.5) +
stat_qq_line() +
stat_qq_point() +
labs(x = "理论分位数", y = "抽样分位数") +
scale_fill_discrete("Bandtype") +
theme_minimal()
# 显示图形
grid.arrange(hist_plot, qq_plot, ncol = 2, top = textGrob(paste0(name,"构造的频数图和密度线，K-S检验下p值为：",ks.test(samples, "pnorm", mean = 0, sd = 1)$p.value), gp = gpar(fontsize = 16, fontface = "bold")))
}这里我们展示四个默认参数（累加数n=12）下生成个500样本正态分布情况，值得说明的是在这一参数下每种方法生成的正态分布对应的p值并不稳定，但总体上都落在K-S检验的置信域内，通过以下右图不同方法检验的置信域绘制，可以看到不是所有的基础分布生成的正态分布都落在不同的置信域内：
# plot_distribution(get_norm_by_unif(m=500),name="平均分布")
# plot_distribution(get_norm_by_exp(m=500),name="指数分布")
# plot_distribution(get_norm_by_gamma(m=500),name="伽马分布")
# plot_distribution(get_norm_by_beta(m=500),name="贝塔分布")
我们可以知道，累加次数n应当是影响生成正态随机数质量的重要参数，根据中心极限定理越大的n能够使累加值越收敛于标准正态分布。在近似抽样方法中，使用不同的基础分布累加时，依赖的n显然应该也是不同的，例如我们以下展示的代码，在遍历1~20的累加次数，循环生成50次、每次500个样本，并将50次得到的K-S检验p值储存到对应的n下，于是就可以画出不同累加次数n对应的p值分布的箱线图。
为什么要检测50次？正如我们前面所说，生成的随机数本身的随机性会影响K-S检验的结果，我们不妨多试验几次，将p值的多次检验结果当做某种意义上度量“精确性”的指标的分布。
我们认为直观上，这个方法绘制出的箱线图应该有这样的性质：
随着n增大，p值箱线图均值和整体应该增高或稳定高于0.05
最初n很小时p值箱线图应该在0.05以下，说明生成的随机数和标准正态分布显著不同
不同的基础分布累加方法对应的p值箱线图在n的递增下应该会具有不同的形态。
# 灵敏度分析函数
sensitivity_analysis <- function(max_n = 20, iterations = 50, m = 500) {
# 初始化存储 p 值的列表
p_values_list <- list()
# 循环 n 从 1 到 max_n
for (n in 1:max_n) {
p_values <- numeric(iterations) # 存储当前 n 下的 p 值
for (i in 1:iterations) {
# 生成样本并进行 K-S 检验
samples <- get_norm_by_unif(n = n, m = m)
p_values[i] <- shapiro.test(samples)$p.value
}
p_values_list[[n]] <- p_values
}
p_values_df <- do.call(rbind, lapply(1:max_n, function(n) {
data.frame(n = n, p_value = p_values_list[[n]])
}))
# 绘制箱线图
ggplot(p_values_df, aes(x = factor(n), y = p_value)) +
geom_boxplot(fill = "#3288bd", alpha = 0.7) +
labs(title = "Sensitivity Analysis of p-values by n",
x = "Number of Additive Steps (n)",
y = "p-value") +
theme_minimal()
}
# 运行灵敏度分析
sensitivity_analysis(max_n = 20, iterations = 50, m = 500)
同理，我们不难对其他四个分布均作出该灵敏度分析，于是可以得到下图：
观察四幅图中不难发现，指数分布构造的随机数在n相对大时仍然有部分会被拒绝，即在95%置信水平下显著与标准正态不同，伽马分布也相对不是很好，贝塔分布相对好一点，均匀分布最佳。通过前面基础分布的频数图绘制的再观察，我们可以发现结论：带有偏态的基础分布会影响累加收敛到正态分布的速度，偏度越高影响越大。
这很好解释了为什么在大多情况下通常是采用均匀分布的累加来构造标准正态随机数的，不过均匀分布的唯一一个限制是，在固定n的情况下，生成随机数的范围限制在\([-n,n]\)之间，这和标准正态分布的范围\((-\infty,\infty)\)不符。
最后我们还意识到一个问题，那就是当我们用某种方法生成随机数已经足够好时，即该方法生成的随机数在K-S检验下的p值基本大于0.05的情况，p值将是一个与“抽样方法”关系不大、单纯因为“生成随机性”而变化的值，并不能很好地进一步刻画某种方法生成随机数的质量。
换句话说，利用假设检验度量“随机数的质量”的方法大多数情况下只能作为一个“是否通过”的保障，在p值大于0.05后就无法获取更多关于该方法好坏的信息了。
在这一部分，我们关心的是，不同基础分布的参数是否会影响生产随机数的质量？
我们以两个分布为例：
1.指数分布
通过该图随机模拟的结果，可以大致认为指数分布\(Exp(\lambda)\)的参数\(\lambda\)对生成的正态分布的质量影响不大。
2.贝塔分布
我们只调节了贝塔分布\(Beta(alpha,beta)\)的\(\beta\)参数，并在右侧一并画出不同参数下对应的基础分布的图像，这样我们就可以观察在不同的均值\(\mu=EX=\frac{a}{a+b}\)下由beta分布作为基础分布生成正态随机数的质量。
在上图左侧子图中可以看到，随着beta越来越大，p值箱线图总体大于0.05所要求的最小累加值n似乎逐渐变大，这意味着相对于固定的alpha值，越大的beta将降低相同累加值下生成随机数的质量。
这是为什么呢？首先beta的增大会导致均值“比率”意义的减少，对于偏度\(\beta_s\)和峰度\(\beta_k\)： \[ \beta_s=\frac{\mathrm{E}(X-\mu)^{3}}{\left[\mathrm{E}(X-\mu)^{2}\right]^{3 / 2}}=\frac{2(b-a) \sqrt{a+b+1}}{(a+b+2) \sqrt{a b}} \\\beta_k=\frac{\mathrm{E}(X-\mu)^{4}}{\left[\mathrm{E}(X-\mu)^{2}\right]^{2}}-3=\frac{6\left[a^{3}-a^{2}(2 b-1)+b^{2}(b+1)-2 a b(b+2)\right]}{a b(a+b+2)(a+b+3)} \\ \] 直接显示求出单调性略微困难，我们不妨直接控制\(a=2\)，绘制出偏度和峰度随着\(b\)的变化图像：
# 定义参数 a
a <- 2
# 定义 b 的范围
b_values <- seq(0.1, 20, by = 0.1)
beta_s <- numeric(length(b_values))
beta_k <- numeric(length(b_values))
for (i in 1:length(b_values)) {
b <- b_values[i]
# 计算偏度
beta_s[i] <- (2 * (b - a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b))
# 计算峰度
beta_k[i] <- (6 * (a^3 - a^2 * (2 * b - 1) + b^2 * (b + 1) - 2 * a * b * (b + 2))) /
(a * b * (a + b + 2) * (a + b + 3)) - 3
}
par(mfrow = c(1, 2))
plot(b_values, beta_s, type = "l", col = "blue",
xlab = "b 值", ylab = "偏度 β_s",
main = "偏度随 b 变化的图像")
plot(b_values, beta_k, type = "l", col = "red",
xlab = "b 值", ylab = "峰度 β_k",
main = "峰度随 beta 变化的图像")
可以看到偏度在b>a=2时逐渐增大偏离正态分布，而峰度在b>2=2时先增后减，在这一趋势下，beta分布产生随机数的质量逐渐变差，也就说，偏度的偏倚情况的影响会略大于峰度，这和我们之前的结论仍然一致。
最后我们来比较这四种方法的时间效率，我们可以通过时间和空间复杂度来进行比较：
library(microbenchmark)
library(ggplot2)
library(pryr)
# 记录运行时间和内存消耗
benchmark_memory <- function(func) {
mem_usage <- mem_change(func()) # 记录内存消耗
return(mem_usage)
}
# 记录四个函数的时间和内存消耗
results <- microbenchmark(
"均匀分布" = { mem_usage_unif <- benchmark_memory(get_norm_by_unif) },
"指数分布" = { mem_usage_exp <- benchmark_memory(get_norm_by_exp) },
"伽马分布" = { mem_usage_gamma <- benchmark_memory(get_norm_by_gamma) },
"贝塔分布" = { mem_usage_beta <- benchmark_memory(get_norm_by_beta) },
times = 100
)
memory_usages <- c(
"均匀分布" = mem_usage_unif,
"指数分布" = mem_usage_exp,
"伽马分布" = mem_usage_gamma,
"贝塔分布" = mem_usage_beta
)
print(memory_usages)均匀分布 指数分布 伽马分布 贝塔分布
736 736 736 736 总体来说内存消耗一致。
library(dplyr)
Attaching package: 'dplyr'The following object is masked from 'package:pryr':
whereThe following object is masked from 'package:gridExtra':
combineThe following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionresults_df <- as.data.frame(results)
# 删除时间超过400毫秒（400000微秒）的异常值
results_df <- results_df %>%filter(time / 1e6 <= 400)
ggplot(results_df, aes(x = expr, y = time / 1e6, fill = expr)) +
geom_boxplot() +
scale_fill_manual(values = c("#356d67", "#4c9568", "#7fb961", "#b0d45d")) +
theme_minimal() +
labs(title = "不同基础分布生成正态随机数的耗时",
x = "分布",
y = "时间 (ms)") +
theme(legend.position = "none")
图中可以发现四种方法耗时相差不大，但最多异常情况的方法指数分布>均匀分布>贝塔分布>伽马分布。
这里我们采用将概率密度函数分解并利用线性近似的Bulter抽样法，初始设置为\(m=40\)个小区间，\(n=500\)个样本，构造函数如下：
get_norm_by_pdf <- function(n = 40,m=500){
eta = NULL
Npdf = NULL
Ncdf = NULL
fx <- function(x) {
return(dnorm(x)) # 使用内置的dnorm函数
}
Fx <- function(x) {
return(pnorm(x)) # 使用内置的pnorm函数
}
x = seq(-6, 6, 0.0001)
Npdf = fx(x)
for (i in 1:length(x)) {
Ncdf[i] = Fx(x[i])
}
xd = NULL
for (i in 1:n) {
# 确定m个小区间的分割点
index = min(which(abs(Ncdf - i/n) < 0.0001))
xd[i] = x[index]
}
xd = c(x[1], xd)
for (i in 1:m) {
r = runif(1, 0, 1)
k = floor(n * r + 1)
u1 = runif(1, 0, 1)
u2 = runif(1, 0, 1)
dk = abs(fx(xd[k+1]) - fx(xd[k])) / (fx(xd[k+1]) + fx(xd[k]))
if (u2 <= dk) {
if (fx(xd[k+1]) > fx(xd[k])) {
eta[i] = xd[k] + (xd[k+1] - xd[k]) * sqrt(u1)
} else {
eta[i] = xd[k+1] - (xd[k+1] - xd[k]) * sqrt(1-u1)
}
} else {
eta[i] = xd[k] + (xd[k+1] - xd[k]) * u1
}
}
return(eta)
}
eta <- get_norm_by_pdf()
# 频率直方图和概率密度函数
hist(eta, 40, xlim=c(-6, 6), freq=FALSE)
x = seq(-6, 6, 0.0001)
lines(x, dnorm(x, 0, 1), lwd=1, col='red')
在以上构造的函数中，主要变动的参数就是区间分割的数量，显然\(n\)越大，分割的区间越多，近似效果越好，我们可以通过以下代码来观察不同\(n\)下的效果，我们遍历5~60，以3为步长的分割数：
max_n = 50
iterations = 20
m = 500
# 初始化存储 p 值的列表
all_n <- floor((max_n-5)/3)+1
p_values_list <- vector("list",all_n)
plots_list <- list()
j = 1
# 循环 n 从 5 到 max_n
for (n in seq(5, max_n, by=3)) {
# print(n)
p_values <- numeric(iterations) # 存储当前 n 下的 p 值
for (i in 1:iterations) {
# 生成样本并进行 K-S 检验
samples <- get_norm_by_pdf(n = n, m = m)
p_values[i] <- ks.test(samples, "pnorm", mean = 0, sd = 1)$p.value
}
# 绘制QQ图
qq_plot <- ggplot(data = data.frame(samples), mapping = aes(sample = samples)) +
geom_qq_band(bandType = "ks", mapping = aes(fill = "KS"), alpha = 0.5) +
geom_qq_band(bandType = "ts", mapping = aes(fill = "TS"), alpha = 0.5) +
geom_qq_band(bandType = "pointwise", mapping = aes(fill = "Normal"), alpha = 0.5) +
geom_qq_band(bandType = "boot", mapping = aes(fill = "Bootstrap"), alpha = 0.5) +
stat_qq_line() +
stat_qq_point() +
labs(x = "理论分位数", y = "抽样分位数") +
scale_fill_discrete("Bandtype") +
theme_minimal()
plots_list[[j]] <- qq_plot
p_values_list[[j]] <- p_values
j = j + 1
}
# 将结果转换为数据框以便绘图
p_values_df <- do.call(rbind, lapply(1:all_n, function(n) {
data.frame(n = n, p_value = p_values_list[[n]])
}))
# 绘制箱线图
ggplot(p_values_df, aes(x = factor(n), y = p_value)) +
geom_boxplot(fill = "#3288bd", alpha = 0.7) +
labs(title = "Sensitivity Analysis of p-values by n",
x = "Number of Segments (n)",
y = "p-value") +
#添加p = 0.05的参考线，并标注数值
geom_hline(yintercept = 0.05, linetype = "dashed", color = "red") +
annotate("text", x = 0, y = 0.1, label = "p = 0.05", color = "red", hjust = 0, vjust = 0) +
theme_minimal()grid.arrange(plots_list[[1]],plots_list[[8]],plots_list[[16]],ncol = 3)
我们在这里特殊说明，由于为了统计方便，这里横坐标的n并不是实际分割区间的数量，实际分割区间数量与这里的n对应的关系为\(n_{true}=5+3\cdot (n-1)\)这一部分我们可以看到，当\(n_{true}\)越大时，生成的随机数质量相对越好。
我们抽取\(n_{true}=5,26,50\)三种分割区间数量下的QQ图，同样可以看到，当分割区间数量越多时，生成的随机数越接近标准正态分布在QQ图上呈现的直线，当分割区间较少时生成的随机数相对弯曲且在部分检验方法的置信域外。
library(microbenchmark)
library(pryr)
# 记录运行时间和内存消耗
benchmark_memory <- function(func,n) {
mem_usage <- mem_change(func(n=n)) # 记录内存消耗
return(mem_usage)
}
# 记录四个函数的时间和内存消耗
results <- microbenchmark(
"n=5" = { mem_usage_5 <- benchmark_memory(get_norm_by_pdf,n = 5) },
"n=10" = { mem_usage_10 <- benchmark_memory(get_norm_by_pdf,n = 10) },
"n=15" = { mem_usage_15 <- benchmark_memory(get_norm_by_pdf,n = 15) },
"n=20" = { mem_usage_20 <- benchmark_memory(get_norm_by_pdf,n = 20) },
"n=25" = { mem_usage_25 <- benchmark_memory(get_norm_by_pdf,n = 25) },
"n=30" = { mem_usage_30 <- benchmark_memory(get_norm_by_pdf,n = 30) },
"n=35" = { mem_usage_35 <- benchmark_memory(get_norm_by_pdf,n = 35) },
"n=40" = { mem_usage_40 <- benchmark_memory(get_norm_by_pdf,n = 40) },
"n=45" = { mem_usage_45 <- benchmark_memory(get_norm_by_pdf,n = 45) },
"n=50" = { mem_usage_50 <- benchmark_memory(get_norm_by_pdf,n = 50) },
"n=55" = { mem_usage_55 <- benchmark_memory(get_norm_by_pdf,n = 55) },
"n=60" = { mem_usage_60 <- benchmark_memory(get_norm_by_pdf,n = 60) },
times = 30
)
memory_usages <- c(
"n=5" = mem_usage_5,
"n=10" = mem_usage_10,
"n=15" = mem_usage_15,
"n=20" = mem_usage_20,
"n=25" = mem_usage_25,
"n=30" = mem_usage_30,
"n=35" = mem_usage_35,
"n=40" = mem_usage_40,
"n=45" = mem_usage_45,
"n=50" = mem_usage_50,
"n=55" = mem_usage_55,
"n=60" = mem_usage_60
)
print(memory_usages) n=5 n=10 n=15 n=20 n=25 n=30 n=35 n=40 n=45 n=50 n=55 n=60
736 736 736 736 736 736 736 736 736 736 736 736 results_df <- as.data.frame(results)
# 删除时间超过400毫秒（400000微秒）的异常值
# results_df <- results_df %>% filter(time / 1e6 <= 400)
ggplot(results_df, aes(x = expr, y = time / 1e6, fill = expr)) +
geom_boxplot() +
scale_fill_brewer(palette = "Set3") + # 使用离散调色板
theme_minimal() +
labs(title = "不同分割区间数目生成正态随机数的耗时",
x = "分割区间数目=5+3(n-1)",
y = "时间 (ms)") +
theme(legend.position = "none")
我们同样对不同分割区间数量下的时间效率进行了比较，可以看到，随着分割区间数量的增加，生成随机数的时间消耗的均值基本呈现增加趋势，这是因为除了区间分割的判断需要\(n\)次循环，还要在每次生成随机数时都有重新进行循环抽样，这一部分的时间复杂度为\(O(m+n)\)。
将该算法和前面的中心极限定理累加的近似方法对比，空间消耗上基本一致均为736，可能因为内存大小主要受最终抽样数量\(n=500\)决定，其中中间过程产生的变量影响不大；时间效率上，后者耗时四个方法基本在400ms以内，几乎为该算法的一半，这是因为中心极限定理的累加方法的时间复杂度为\(O(1)\)，即不随着累加次数的增加而增加，主要的时间消耗只在根据设置的累加次数要一次性生成多少个基础分布的随机数。
get_norm_by_cdf <- function(n = 500, m = 500) {
# 生成n个标准正态分布的观测数据
observed_data <- rnorm(n, mean = 0, sd = 1)
# 对观测数据进行排序
sorted_data <- sort(observed_data)
# 生成m个随机数用于抽样
eta <- numeric(m)
for (i in 1:m) {
# 生成U1 ~ U(0, 1)
U1 <- runif(1, min = 0, max = 1)
# 找到对应的区间 (x(k), x(k+1))
k <- floor(U1 * (n - 1)) + 1
# 生成U2 ~ U(0, 1) 并计算eta
U2 <- runif(1, min = 0, max = 1)
eta[i] <- sorted_data[k] + U2 * (sorted_data[k+1] - sorted_data[k])
}
return(eta)
}
# 调用函数，生成抽样数据
eta <- get_norm_by_cdf()
# 绘制频率直方图和理论的正态分布曲线
hist(eta, 40, xlim = c(-6, 6), freq = FALSE, main = "抽样数据的频率直方图")
x <- seq(-6, 6, length = 100)
lines(x, dnorm(x, mean = 0, sd = 1), lwd = 1, col = 'red')
max_n = 1000
iterations = 20
m = 500
# 初始化存储 p 值的列表
all_n <- floor((max_n-5)/100)+1
p_values_list <- vector("list",all_n)
plots_list <- list()
j = 1
# 循环 n 从 5 到 max_n
for (n in seq(100, max_n, by=100)) {
# print(n)
p_values <- numeric(iterations) # 存储当前 n 下的 p 值
for (i in 1:iterations) {
# 生成样本并进行 K-S 检验
samples <- get_norm_by_cdf(n = n, m = m)
p_values[i] <- ks.test(samples, "pnorm", mean = 0, sd = 1)$p.value
}
# 绘制QQ图
qq_plot <- ggplot(data = data.frame(samples), mapping = aes(sample = samples)) +
geom_qq_band(bandType = "ks", mapping = aes(fill = "KS"), alpha = 0.5) +
geom_qq_band(bandType = "ts", mapping = aes(fill = "TS"), alpha = 0.5) +
geom_qq_band(bandType = "pointwise", mapping = aes(fill = "Normal"), alpha = 0.5) +
geom_qq_band(bandType = "boot", mapping = aes(fill = "Bootstrap"), alpha = 0.5) +
stat_qq_line() +
stat_qq_point() +
labs(x = "理论分位数", y = "抽样分位数") +
scale_fill_discrete("Bandtype") +
theme_minimal()
plots_list[[j]] <- qq_plot
p_values_list[[j]] <- p_values
j = j + 1
}
# 将结果转换为数据框以便绘图
p_values_df <- do.call(rbind, lapply(1:all_n, function(n) {
data.frame(n = n, p_value = p_values_list[[n]])
}))
# 绘制箱线图
ggplot(p_values_df, aes(x = factor(n), y = p_value)) +
geom_boxplot(fill = "#3288bd", alpha = 0.7) +
labs(title = "Sensitivity Analysis of p-values by n",
x = "Number of Segments (n*100)",
y = "p-value") +
#添加p = 0.05的参考线，并标注数值
geom_hline(yintercept = 0.05, linetype = "dashed", color = "red") +
annotate("text", x = 0, y = 0.1, label = "p = 0.05", color = "red", hjust = 0, vjust = 0) +
theme_minimal()
基本可以看到，随着观测数\(n\)的增加，生成的随机数质量越好，但是在这个样本观测值范围下，生成随机数的稳定性不如根据上述参数设置的密度函数抽样方法，在观测数\(n=1000\)时，仍有小部分抽样数据的p值小于0.05。
选取\(n=100,500,1000\)画出Q-Q图，可以观察到100时仍然不太稳定，500后已经基本都落在不同检验方法的置信区域内。
grid.arrange(plots_list[[1]],plots_list[[5]],plots_list[[10]],ncol = 3)library(microbenchmark)
library(pryr)
# 记录运行时间和内存消耗
benchmark_memory <- function(func,n) {
mem_usage <- mem_change(func(n=n)) # 记录内存消耗
return(mem_usage)
}
# 记录四个函数的时间和内存消耗
results <- microbenchmark(
"n=100" = { mem_usage_100 <- benchmark_memory(get_norm_by_cdf,n = 100) },
"n=200" = { mem_usage_200 <- benchmark_memory(get_norm_by_cdf,n = 200) },
"n=300" = { mem_usage_300 <- benchmark_memory(get_norm_by_cdf,n = 300) },
"n=400" = { mem_usage_400 <- benchmark_memory(get_norm_by_cdf,n = 400) },
"n=500" = { mem_usage_500 <- benchmark_memory(get_norm_by_cdf,n = 500) },
"n=600" = { mem_usage_600 <- benchmark_memory(get_norm_by_cdf,n = 600) },
"n=700" = { mem_usage_700 <- benchmark_memory(get_norm_by_cdf,n = 700) },
"n=800" = { mem_usage_800 <- benchmark_memory(get_norm_by_cdf,n = 800) },
"n=900" = { mem_usage_900 <- benchmark_memory(get_norm_by_cdf,n = 900) },
"n=1000" = { mem_usage_1000 <- benchmark_memory(get_norm_by_cdf,n = 1000) },
times = 100
)
memory_usages <- c(
"n=100" = mem_usage_100,
"n=200" = mem_usage_200,
"n=300" = mem_usage_300,
"n=400" = mem_usage_400,
"n=500" = mem_usage_500,
"n=600" = mem_usage_600,
"n=700" = mem_usage_700,
"n=800" = mem_usage_800,
"n=900" = mem_usage_900,
"n=1000" = mem_usage_1000
)
print(memory_usages) n=100 n=200 n=300 n=400 n=500 n=600 n=700 n=800 n=900 n=1000
736 736 736 736 736 736 736 736 736 736 results_df <- as.data.frame(results)
# 删除时间超过400毫秒（400000微秒）的异常值
# results_df <- results_df %>% filter(time / 1e6 <= 400)
ggplot(results_df, aes(x = expr, y = time / 1e6, fill = expr)) +
geom_boxplot() +
scale_fill_brewer(palette = "Set3") + # 使用离散调色板
theme_minimal() +
labs(title = "不同观测数据数目生成正态随机数的耗时",
x = "观测数据数目=100n",
y = "时间 (ms)") +
theme(legend.position = "none")
这里的内存消耗仍然是736和前面的解释一致；时间效率上，该算法和前面通过密度函数抽样的方法不同，因为尽管算法内的循环次数对应的时间复杂度为仍为\(O(m)\)，但观测数据量\(n\)的改变只对生成观测数据的一步产生了几乎可以忽略不计的影响，所以任意观测数据量下的时间消耗基本一致，且均在300ms（均值目测150ms）以内，比中心极限定理的累加方法略高。
在本次实验中，我们通过中心极限定理的累加方法、密度函数近似抽样方法和经验分布近似抽样方法生成了500个正态分布的随机数，并通过K-S检验和QQ图对生成的随机数质量进行了评估，最后对不同参数的影响和不同方法的效率进行了灵敏度分析。
综合来看，利用中心极限定理的累加方法生成的随机数质量相对较好，且时间效率较高，但是对于不同基础分布的种类和对应的参数变化较为敏感，通常情况下和正态分布偏度越接近的基础分布效果更好，即能在同等累加次数下更快的通过正态性检验；
密度函数近似抽样方法生成的随机数质量次之，但要求较高的区间分割数目（40-60），耗费时间较长且小幅随分割区间数目增大而增大；
经验分布近似抽样方法生成的随机数质量在观测样本足够多时效果不错，其优势主要在于相对密度函数抽样算法流程较简单，时间快（略慢于中心极限定理累加方法）且几乎不受观测样本量的影响。
简单对比三种方法的优劣表格罗列如下：
| 方法名称 | 生成质量 | 时间效率 | 空间效率 |
|---|---|---|---|
| 中心极限定理累加方法 | 较好 | 高 | 一般 |
| 密度函数近似抽样方法 | 一般 | 低 | 一般 |
| 经验分布近似抽样方法 | 较好 | 中 | 一般 |
题目：
Answer：
要求生成三维随机向量 \((x_1, x_2, x_3)\)，其密度函数为：
\[ f(x_1, x_2, x_3) = \frac{4}{3\pi} \quad \text{对于} \quad x_1^2 + x_2^2 + x_3^2 < 1 \]
这是一个在单位球面内的均匀分布问题，注意到\(f_0 = \sup_{x_1,x_2,x_3} f(x_1,x_2,x_3) = \frac{4}{3\pi}\)，且\(-1<x_i<1,i=1,2,3\)
step1：生成三个在\([0,1]\)上的均匀分布的随机数\(U_0,U_1,U_2,U_3\)
step2：计算\(X_1 = 2U_1-1,X_2 = 2U_2-1,X_3 = 2U_3-1\)，则\(X_1 \sim U(-1,1),X_2 \sim U(-1,1),X_3 \sim U(-1,1)\)
step3：如果\(X_1^2 + X_2^2 + X_3^2 < 1\)，则输出\((X_1,X_2,X_3)\)，否则返回step1
伪代码：
while(TRUE):
U_i = GenerateUniformRandomNumber(0, 1),i = 0,1,2,3
X_j = 2 * U_j - 1,j=1,2,3
if (X_1^2 + X_2^2 + X_3^2) < 1 :
return (X_1, X_2, X_3)
else:
continue题目：
设坛子中有 \(n\)个不同颜色的球，共计 \(r\)种颜色，其中第 \(i\)种颜色的球有 \(n_i\) 个，\(n=\sum_{i=1}^r n_i\)。从坛子中随机无放回地抽取 \(m\)个球，设随机变量 \(X_i\) 表示取出的第 i种颜色的球的个数，试产生随机向量 \(\left(X_1, X_2,\cdots, X_r\right)^T\) 的随机数。
Answer：
该题显然适合用条件分布抽样法，\(X_1\)的分布为超几何分布，即：
\[ P(X_1 = k) = \frac{C_{n_1}^k C_{n-n_1}^{m-k}}{C_n^m},k=0,1,\cdots,min(n_1,m) \] 给定\(X_1 = x_1\)后，\(X_2\)的分布为超几何分布，即：
\[ P(X_2 = k|X_1 = x_1) = \frac{C_{n_2}^k C_{n-n_1-n_2}^{m-x_1-k}}{C_{n-n_1}^{m-x_1}},k=0,1,\cdots,min(n_2,m-x_1) \]
以此类推，给定\((X_1,X_2,\cdots,X_{i-1})=(x_1,x_2,\cdots,x_{i-1})\)后，\(X_i\)的分布为超几何分布，即：
\[ P(X_i = k|(X_1,X_2,\cdots,X_{i-1})) = \frac{C_{n_i}^k C_{n-\sum_{j=1}^{i-1}n_j}^{m-\sum_{j=1}^{i-1}x_j-k}}{C_{n-\sum_{j=1}^{i-1}n_j}^{m-\sum_{j=1}^{i-1}x_j}},k=0,1,\cdots,min(n_i,m-\sum_{j=1}^{i-1}x_j) \] 于是产生随机数可按照以下流程：
step1：产生随机数\(X_1\)的超几何分布
step2：给定\(X_1,X_2,\cdots,X_{i-1}\)后，产生随机数\(X_i\)的超几何分布
伪代码：
初始化：
n_i = N_i ,i = 1,2,...,r
sum_x = 0
sum_ni = n_1
抽样：
for i = 1 to r:
n_max = min(n_i,m-sum_x)
prob_distribution = dhyper(0:n_max,n_i,n-sum_ni,m-sum_x)
X_i = sample(0:n_max, size = 1, prob = prob_distribution)
sum_x = sum_x + X_i
sum_ni = sum_ni + n_i
X = (X_1,X_2,...,X_r)