R 中的光谱测试
Spectral Test in R
我想更改我的代码，将 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)
你所拥有的称为一般形式的线性同余生成器
Xn+1 = (aXn+c) mod m
目标是选择这样的 a、c 和 m，使生成的序列尽可能类似于随机序列。没有选择 a、c 和 m 的最佳方法，但有些事情我们可以轻松做到。
你已经选择了 m = 10000。那么我们可以使用一个众所周知的定理（参见，例如，这个 post 的结尾）如何选择这样的 a 和 c，使生成的数字开始仅在 m 步后重复自己。
条件1：c和m需互质。我们有 m = 10000 = 2454。同时，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。
我想更改我的代码，将 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)
你所拥有的称为一般形式的线性同余生成器
Xn+1 = (aXn+c) mod m
目标是选择这样的 a、c 和 m，使生成的序列尽可能类似于随机序列。没有选择 a、c 和 m 的最佳方法，但有些事情我们可以轻松做到。
你已经选择了 m = 10000。那么我们可以使用一个众所周知的定理（参见，例如，这个 post 的结尾）如何选择这样的 a 和 c，使生成的数字开始仅在 m 步后重复自己。
条件1：c和m需互质。我们有 m = 10000 = 2454。同时，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。