r - R中的光谱测试

Question

我想以将 10000 个点放在 [0,1]^2 图中的方式更改我的代码。我尝试将 256 更改为 10000，但它会产生奇怪的位置。我应该改变因子 137 和 187 但不知道如何改变它。有人知道背后的逻辑吗？

工作样本：

nSim = 256
X=rep(0,nSim)
for (i in 2:nSim){
    X[i] = (137*X[i-1]+187)%%256 
}
plot(X[-1],X[-nSim],col="blue",type="p",pch="x",lwd=2)

我的代码：

nSim = 10000
X=rep(0,nSim)
for (i in 2:nSim){
  X[i] = ((137*X[i-1]+187)%%nSim)
}
plot(X[-1]/nSim,X[-nSim]/nSim,col="blue", type="p",pch=20,lwd=2)

score 3 · Accepted Answer

你所拥有的称为线性同余生成器，采用一般形式

X _n+1 = (aX _n +c) mod m

目标是选择这样的 a、c 和 m，以使生成的序列尽可能地类似于随机序列。没有选择 a、c 和 m 的最佳方法，但我们可以轻松做到。

您已经选择了 m = 10000。然后我们可以使用一个众所周知的定理（例如，参见本文的结尾）如何选择这样的 a 和 c，以使生成的数字仅在 m 步后才开始重复。

条件 1：c 和 m 需要互质。我们有 m = 10000 = 2 ⁴ 5 ⁴。同时，187 既不能被 2 也不能被 4 整除，所以我们很好。

条件2：如果4整除m，那么4也需要整除a-1。在我们的例子中，137-1 = 136 可以被 4 整除：136/4 = 34，所以我们在这里也很好。

条件3：如果任意素数p整除m，那么p也需要整除a-1。我们已经在上一步中检查了 p=2，所以剩下 p=5。但是 5 不能除以 137-1 = 136！确实，因此我们得到

length(unique(X))
# [1] 2000

也就是说，在长度为 10000 的序列中，我们只有 2000 个唯一数字，这意味着相同的数字重复了五次。所以，那么我们需要这样一个 a-1 可以被 4 和 5 整除。这提供了很多选择！例如，我们可以选择 a = 4*5*6 + 1 = 121。

因此，使用 a = 121、c = 187 和 m = 10000 给出

length(unique(X))
# [1] 10000

和一个情节

它仍然有些常规，但绝对比以前的要好。您可以继续尝试满足这三个条件的不同 a 和 c。

1 回答 1