Неужели LCG провалил тест Колмогорова-Смирнова так сильно, как подсказывает мой код? - PullRequest
Новая
       23

Неужели LCG провалил тест Колмогорова-Смирнова так сильно, как подсказывает мой код?

1 голос
/ 12 января 2020

Я использую следующий код Python, чтобы проиллюстрировать генерацию случайных переменных для студентов: 

import numpy as np
import scipy.stats as stats

def lcg(n, x0, M=2**32, a=1103515245, c=12345):
    result = np.zeros(n)
    for i in range(n):
        result[i] = (a*x0 + c) % M
        x0 = result[i]

    return np.array([x/M for x in result])

x = lcg(10**6, 3)
print(stats.kstest(x, 'uniform'))

Параметры по умолчанию - те, что используются glib c, согласно Википедии. Последняя строка кода печатает 

KstestResult(statistic=0.043427751892089805, pvalue=0.0)

Значение 0,0 указывает, что наблюдение в принципе никогда не будет происходить, если бы элементы x были действительно распределены в соответствии с равномерным распределением. Мой вопрос: есть ли ошибка в моем коде, или LCG с заданными параметрами не проходит тест Колмогорова-Смирнова с 10**6 репликами?

Ответы [ 2 ]

2 голосов
/ 12 января 2020

Есть проблема с вашим кодом, он делает равномерное распределение как

enter image description here

Я немного изменил вашу реализацию LCG, и все хорошо сейчас (Python 3,7, Анаконда, Win10 x64) 

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

def lcg(n, x0, M=2**32, a=1103515245, c=12345):
    result = np.zeros(n)
    for i in range(n):
        x0 = (a*x0 + c) % M
        result[i] = x0

    return np.array([x/float(M) for x in result])

#x = np.random.uniform(0.0, 1.0, 1000000)
x = lcg(1000000, 3)
print(stats.kstest(x, 'uniform'))

count, bins, ignored = plt.hist(x, 15, density=True)
plt.plot(bins, np.ones_like(bins), linewidth=2, color='r')
plt.show()

, которая печатает 

KstestResult(statistic=0.0007238884545415214, pvalue=0.6711878724246786)

и графики

enter image description here

ОБНОВЛЕНИЕ

, как указывало @ p js, лучше делить на float (M) прямо в l oop, нет необходимости во втором проходе по всему массиву 

def lcg(n, x0, M=2**32, a=1103515245, c=12345):
    result = np.empty(n)
    for i in range(n):
        x0 = (a*x0 + c) % M
        result[i] = x0 / float(M)

    return result
1 голос
/ 13 января 2020

В дополнение к ответу Северина причина, по которой мой код не работал должным образом, заключается в том, что result был массивом чисел с плавающей запятой. Мы видим разницу между двумя реализациями уже на второй итерации. После первой итерации x0 = 3310558080. 

In [9]: x0 = 3310558080

In [10]: float_x0 = float(x0)

In [11]: (a*x0 + c) % M
Out[11]: 465823161

In [12]: (a*float_x0 + c) % M
Out[12]: 465823232.0

In [13]: a*x0
Out[13]: 3653251310737929600

In [14]: a*float_x0
Out[14]: 3.6532513107379297e+18

Таким образом, проблема была связана с использованием чисел с плавающей запятой.

...