来自 [0, 1] 上的均匀分布的随机样本答案

【问题标题】：Random samples from a uniform distribution over [0, 1]来自 [0, 1] 上的均匀分布的随机样本
【发布时间】：2016-01-21 06:55:19
【问题描述】：

我们都知道并喜欢numpy 函数random.rand，它“创建[s] 一个给定形状的数组并使用来自[0, 1) 上的均匀分布的随机样本传播[s]”：

>>> import numpy as np
>>> np.random.rand(3,2)
array([[ 0.14022471,  0.96360618],  #random
       [ 0.37601032,  0.25528411],  #random
       [ 0.49313049,  0.94909878]]) #random

如果我想要来自[0, 1] 上的均匀分布的随机样本，我有哪些选择？ 注意细微差别： numpy 函数不包括 1；我想要的功能包括 1.

【问题讨论】：

在 float64 中，如果要包含完美的 1，则恰好有 2^52 分之一的机会发生。我的问题是，这真的重要吗？
另外，如果您要复制论文，您可能想弄清楚它是使用 Mersenne Twister（本机 python 实现）还是实际的真随机数生成器。
哈哈，是的。尝试复制实验结果始终是个问题。
我之前曾尝试回答您的问题（见下文），但与其他人一样，我强烈认为您不应该关心这个非常微小的差异。我还认为，如果这种差异对您很重要，您是否知道非规范化浮点数 en.wikipedia.org/wiki/Denormal_number ？简而言之，编码浮点数的密度在 [0,1] 中的各处并不均匀。我不确切知道 Numpy 是如何生成随机数的，但这样的事实也会让你感到不安。
难点：你为什么特别关心1.0被排除在外？在[0.0, 1.0] 范围内大约有 2^62 个不同的双精度浮点数。 random.rand 已经只能生成其中的一小部分——大约 500 分之一。（random.rand 有 2^53 个可能的输出。）再缺少一个可能的输出值确实不会产生太大的影响。 IOW，只需使用[0, 1) 范围；它与人工创建的[0, 1] 随机范围在统计上无法区分。

标签： python numpy random probability

【解决方案1】：

如果我没记错的话，你可以用类似的方法构建小于 1 的最大数：

import struct
m = struct.unpack("d",struct.pack("Q",eval("0b00" + "1"*9 + "0" + "1"*52)))[0]

这似乎是：

m = struct.unpack("d",'\xff\xff\xff\xff\xff\xff\xef?')[0]

（根据 IEEE 754 标准）。看看：

from decimal import Decimal
Decimal(m)

哪个打印：

Decimal('0.99999999999999988897769753748434595763683319091796875')

然后，只需将您的np.random.rand 除以该值，您就应该得到您想要的。

【讨论】：

直观地说，除法后的舍入会导致其他浮点数被排除在可能的范围之外，不是吗？尽管如此，技术上还是非常正确的。
我想那我还得确认numpy.random.rand实际上包括这个最大可能的小于一的数字？
获得小于 1 的最大浮点数的更简单方法是numpy.nextafter(1.0, 0)。
是的，这看起来肯定会让事情变得更糟。
@Baruchel：如果你假设 IEEE 754，你不妨直接使用 m = 1 - 2**-53 并避免使用 struct。

【解决方案2】：

我之前回答过您的问题（请参阅我的其他回答）；我也发表了评论。正如我在评论中所说，我真的不相信你应该关心这个小区别。但也许你可以像那样实现你想要的。

得到一个随机INTEGER，对应于[1,2]中的一个IEEE 754数的位序列而不是[0,1]因为非规范化的数字；将这些位转换为浮点 IEEE 754 数字；然后减去 1。这将起作用，因为 IEEE 754 中 [1,2] 中的所有数字都是连续的。

# sequence of bit for 1.0
>>> eval("0b00" + "1"*10 + "0"*51 + "0")
4607182418800017408

# sequence of bit for first number after 2.0
>>> eval("0b01" + "0"*10 + "0"*51 + "1")
4611686018427387905

>>> struct.unpack("d",struct.pack("Q",np.random.randint(4607182418800017408,4611686018427387905)))[0]-1
0.9349236864324189

你应该能够得到 0 和 1。但是我不确定两个大整数值之间的整数是否完全一致......

更简洁的功能可能是：

from struct import pack, unpack
import numpy as np

def myrand():
  a = unpack("Q",pack("d",1.0))[0]               # 1.0
  b = unpack("Q",pack("d",np.nextafter(2,3)))[0] # 2.000000000000000444
  return unpack("d",pack("Q",np.random.randint(a,b)))[0]-1

但同样，这更像是一个概念证明。

【讨论】：