Спектральный тест в 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 = 10000 = 2 ⁴ 5 ⁴ .Между тем, 187 не делится ни на 2, ни на 4, поэтому мы хороши.

Условие 2 : если 4 делит m, то 4 также должно делить a-1.В нашем случае 137-1 = 136 делится на 4: 136/4 = 34, поэтому и нам здесь хорошо.

Условие 3 : если любое простое число p делит m, то pтакже необходимо разделить а-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, которые удовлетворяют трем условиям.