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)

Answer 1

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

X_n+1 = (aX_n+c) mod m

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

你已经选择了 m = 10000。那么我们可以使用一个众所周知的定理（例如，参见 post 的结尾）如何选择这样的 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。