LCG 在 Kolmogorov-Smirnov 测试中是不是像我的代码所暗示的那样严重失败？

Posted 2023-03-12

【中文标题】LCG 在 Kolmogorov-Smirnov 测试中是不是像我的代码所暗示的那样严重失败？

【英文标题】：Does the LCG fail the Kolmogorov-Smirnov test as badly as my code suggests?LCG 在 Kolmogorov-Smirnov 测试中是否像我的代码所暗示的那样严重失败？

【发布时间】：2020-04-29 09:45:23

【问题描述】：

我用下面的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'))

根据 Wikipedia，默认参数是 glibc 使用的参数。打印代码的最后一行

KstestResult(statistic=0.043427751892089805, pvalue=0.0)

pvalue 0.0 表示如果x 的元素真正按照均匀分布进行分布，则观察基本上不会发生。 我的问题是：我的代码中是否存在错误，或者具有给定参数的 LCG 是否通过 10**6 副本的 Kolmogorov-Smirnov 测试？

【问题讨论】：

你使用什么版本的 Python？ Python 2 和 Python 3 的除法不同

【参考方案1】：

你的代码有问题，它使分布像

我已经稍微改变了你的 LCG 实现，现在一切都很好（Python 3.7、Anaconda、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)

和情节

更新

正如@pjs 指出的那样，你最好在循环中除以 float(M) ，不需要 第二遍遍历整个数组

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

【参考方案2】：

为了补充 Severin 的回答，我的代码无法正常工作的原因是 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

所以问题与使用浮点数有关。

【讨论】：

您的代码不起作用的原因是因为您使用了x0 = result[i] 而不是result[i] = x0 并且忽略了除以float(M)。

@pjs 他除以 float(M) 但最后效率很低

用更高效的实现更新了我的答案

根据我的实验，除以float(M) 或M 无关紧要。

@Rastapopoulos 在 python 3 中为真，在 python 2 中不为真。当 Severin 询问时，您没有指定使用哪个，float(M) 在这两种情况下都有效。