LCG 是否像我的代码所建议的那样严重未能通过 Kolmogorov-Smirnov 测试?
Does the LCG fail the Kolmogorov-Smirnov test as badly as my code suggests?
我用下面的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'))
根据维基百科,默认参数是 glibc 使用的参数。代码的最后一行打印
KstestResult(statistic=0.043427751892089805, pvalue=0.0)
p值为0.0表示如果x
的元素真正服从均匀分布,则观察基本不会发生。
我的问题是:我的代码中是否存在错误,或者具有给定参数的 LCG 是否未通过 10**6
个副本的 Kolmogorov-Smirnov 测试?
你的代码有问题,它像
一样均匀分布
我稍微更改了您的 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
为了补充 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
所以问题与浮点数的使用有关。
我用下面的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'))
根据维基百科,默认参数是 glibc 使用的参数。代码的最后一行打印
KstestResult(statistic=0.043427751892089805, pvalue=0.0)
p值为0.0表示如果x
的元素真正服从均匀分布,则观察基本不会发生。
我的问题是:我的代码中是否存在错误,或者具有给定参数的 LCG 是否未通过 10**6
个副本的 Kolmogorov-Smirnov 测试?
你的代码有问题,它像
一样均匀分布我稍微更改了您的 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
为了补充 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
所以问题与浮点数的使用有关。