在 R 中生成大小为 n 的完美正态分布样本

Question

我想生成一个 mean = 0、sd = 1 和 size n = 100 的样本，其分布尽可能正态。单独使用 rnorm returns 很多可变性。

我发现的唯一方法是对多个 rnorms 进行平均。

rowMeans(replicate(10000, sort(rnorm(100, 0, 1))))

这个 returns 一个相当令人满意的结果，但我不确定这是最有效的方法。

---

编辑：

我不希望 mean 和 sd 严格等于 0 和 1，而是 "look" 的分布像正态分布（绘制密度曲线时）。

似乎 qnorm 方法比 "average" 方法效果更差：

# qnorm method
x <- qnorm(seq(.00001, .99999, length.out = 100), mean=0, sd=1)
plot(density(x))

# average method
x <- rowMeans(replicate(10000, sort(rnorm(100, mean=0, sd=1))))
plot(density(x))

如果确定性解决方案以更有效的方式返回接近平均方法的结果，我会很高兴。

---

编辑 2：可能的解决方案

根据答案，以下似乎有效，相对于 n 调整边界：

x <- qnorm(seq(1/n, 1-1/n, length.out = n), mean=0, sd=1)

下面是不同 n 值的 qnorm 和 average 方法的比较：

par(mfrow=c(6,2))
for(n in c(10, 20, 100, 500, 1000, 9876)){
  x <- qnorm(seq(1/n, 1-1/n, length.out = n), mean=0, sd=1)
  plot(density(x), col="blue", lwd=2)

  x <- rowMeans(replicate(10000, sort(rnorm(n, mean=0, sd=1))))
  plot(density(x), col="red", lwd=2)
}

Answer 1

如果您想要一个确定性的解决方案，这应该可行

qnorm(seq(0.01, 0.99, length.out = 100))

注意 qnorm(0) 给出 $-\infty$ 而 qnorm(1) 是 $\infty$，所以你需要找到一些合理的界限。

对于 n=100，界限 0.01 和 0.99 似乎效果最好。如果您希望确定性解决方案的边界更远，则需要增加 n.

Answer 2

您是否要创建 100 个具有近似正态分布且均值恰好为零且 sd 恰好为一的数字？这样做：

粗略开始：

> X = rnorm(100)

转移他们：

> X = X-mean(X)

扩展它们：

> X = X/sd(X)

检查一下：

> mean(X)
[1] -7.223497e-18

差不多了

> sd(X)
[1] 1

加油。

这与 scale 函数的作用相同：

> X = rnorm(100)
> mean(X)
[1] -0.007667039
> sd(X)
[1] 0.9336842
> sx = scale(X)
> mean(sx)
[1] 1.437056e-17
> sd(sx)
[1] 1

Answer 3

低差异序列？ halton、faure、sobol、hammersley：示例：

library(randtoolbox)

sequence <-sobol(n=100, dim = 1, init = TRUE, scrambling = 0, seed = 4711, normal = FALSE)
mean(sequence)
[1] 0.4982031
sd(sequence)
[1] 0.2860574

#trial with prng
set.seed(1) 
sequence2 <- runif(100)
mean(sequence2)
[1] 0.5178471
sd(sequence2)
[1] 0.2675848

具有相同点数的低差异序列优于伪随机生成器，请记住，对于均匀随机样本，真实平均值为 0.5，sd 为 0.2886751 (sqrt(1/12))，请查看数字.

(mean(sequence) - 0.5)/0.5   #  -0.0008984375
(mean(sequence2) - 0.5)/0.5  #  -0.008923532
(sd(sequence) - sqrt(1/12))*sqrt(12)
[1] -0.009067992
(sd(sequence2) - sqrt(1/12))*sqrt(12)
[1] -0.07305918

~10倍好，不信换个种子试试

ks.test(sequence,"runif")

    One-sample Kolmogorov-Smirnov test

data:  sequence
D = 0.96268, p-value < 2.2e-16
alternative hypothesis: two-sided

> ks.test(sequence2,"runif")

    One-sample Kolmogorov-Smirnov test

data:  sequence2
D = 0.93956, p-value < 2.2e-16
alternative hypothesis: two-sided

现在进行一些平衡：

    sequence <- c(sequence, 1.0 - sequence)  #balancing the mean = use antithetics
    #or if you want (sequence <- sequence - mean(sequence))
    normal_sample <- qnorm(sequence)
    normal_sample <- normal_sample/sd(normal_sample)
    plot(normal_sample)

Answer 4

您可以使用 bayetestR 包：

library(bayestestR)
x <-  rnorm_perfect(n = 100, mean = 0, sd = 1)
plot(density(x))