包含范围内的随机浮点双精度

Posted 2023-02-25

【中文标题】包含范围内的随机浮点双精度

【英文标题】：Random floating point double in Inclusive Range

【发布时间】：2012-04-01 05:41:08

【问题描述】：

我们可以使用下面列出的函数轻松获得所需范围[X,Y) 内的随机浮点数，因为Math.random()（和大多数伪随机数生成器，AFAIK）在[0,1):

function randomInRange(min, max) 
  return Math.random() * (max-min) + min;

// Notice that we can get "min" exactly but never "max".

我们如何才能在两个边界（即[X,Y]）包含的所需范围内获得随机数？

我想我们可以通过“滚动”IEE-754 floating point double precision 的位将最大可能值精确地设置为 1.0 来从Math.random()（或等效项）“增加”我们的值，但这似乎很难做到正确，特别是在不适合位操作的语言中。有没有更简单的方法？

（顺便说一句，为什么随机数生成器会在 [0,1) 而不是 [0,1] 中生成数字？）

[编辑]请注意，我对此没有需要，并且我完全意识到这种区别是迂腐的。只是好奇并希望得到一些有趣的答案。如果这个问题不合适，请随时投票结束。

【问题讨论】：

您能解释一下 max inclusive 将如何产生任何显着差异吗？这些是浮点数。 0.499999999 应该与 0.5 相差不大。

@KendallFrey：不，我的问题没有实际应用。我只是好奇。如果我只是傻，请随意投票关闭。

见http://***.com/questions/363681/java-generating-random-number-in-a-range

回答你的“旁白”，很可能是因为它们是从整数 PRNG 生成的（将在 0 -> 2^n-1 范围内），然后缩放为浮点数。

【参考方案1】：

考虑到 0 到 1 之间的“非常大”数量的值，这真的很重要吗？ 实际上达到 1 的可能性很小，因此不太可能对您正在做的任何事情产生重大影响。

【讨论】：

不，我的问题没有实际应用。我只是好奇=)

【参考方案2】：

这个问题类似于问，1.0 之前的浮点数是多少？ 有这样一个浮点数，但它是 2^24 分之一（对于 IEEE float ) 或 2^53 分之一（double）。

在实践中差异可以忽略不计。

【讨论】：

@OliCharlesworth：扫荡的是哪一个？

“在实践中差异可以忽略不计”。人们可能会编造一个差异很重要的场景。

@OliCharlesworth 实际上，您可能做不到。如果您对1 和下一个最小浮点数之间的差异很敏感，那么您的结果将被浮点数精度问题所覆盖，并且在您使用某些东西之前不会代表任何有意义的东西具有更高的精度（此时您不再对 1 和下一个最小浮点数之间的差异敏感）。

其实更像是 1/2**53 ;)

@tc：是的，你说的很对。我写得比我想象的要快 2^56，忘记了指数是 11 位，而不是 IEEE-754 shortreal 的 8 位。我已经修复了它，尽管这对我的观点没有影响。

【参考方案3】：

在什么情况下您需要一个浮点值来包含上限？对于我理解的整数，但对于浮点数，包含和排除之间的区别就像 1.0e-32。

【参考方案4】：

这样想。如果您想象浮点数具有任意精度，那么准确得到min 的机会为零。获得max 的机会也是如此。我会让你自己得出结论。

这个'问题'相当于在0到1之间的实线上得到一个随机点。没有'包容'和'排斥'。

【讨论】：

而且得到任何特定值的机会也为零，所以总结起来你不会得到任何结果？

嗯，有（不可数）无限个实数，所以0 * infinity = undefined。实际上，结果是 1。这是因为上式中的 0 是由1(sum of probabilities) / infinity(number of possibilitis) = 0(probability of any one number) 确定的。使用该逻辑，您可以简化0 * inf = 1 / inf * inf = 1。

这个论点的问题是浮点数不具有任意精度！在任何真实世界的实验中，您都可以通过统计方式证明差异。

在任何实际实验中，您都需要花费太多 CPU 年才能得出大量随机数样本。

@KendallFrey：只有 2^32 个不同的单精度浮点数（其中相当一部分不在 [0,1) 范围内）； 20 亿次迭代并没有那么那么长。 （但对于双精度，你可能是安全的。）

【参考方案5】：

我相信有更好的决定，但这个应该可行:)

function randomInRange(min, max) 
  return Math.random() < 0.5 ? ((1-Math.random()) * (max-min) + min) : (Math.random() * (max-min) + min);

【讨论】：

+1 好！除了 0 和 1 的机会都是其余数字的一半。 :)

+1 表示对称，但请注意隐含假设 1-Math.random() 不会丢失精度。

【参考方案6】：

只需选择稍大的半开区间，以便您选择的闭区间是一个子集。然后，继续生成随机变量，直到它落在所述闭区间内。

示例：如果您想要 [3,8] 中的统一随机变量，则在 [3,9) 中重复生成统一随机变量，直到它恰好落在 [3,8] 中。

function randomInRangeInclusive(min,max) 
 var ret;
 for (;;) 
  ret = min + ( Math.random() * (max-min) * 1.1 );
  if ( ret <= max )  break; 
 
 return ret;

注意：生成半开 R.V. 的次数。是随机的并且可能是无限的，但是您可以根据需要使预期的调用次数接近 1，而且我认为不存在不会无限次调用的解决方案。

【讨论】：

它在数学上工作，但在浮点世界中，不能保证它会返回max。

【参考方案7】：

首先，您的代码有问题：尝试randomInRange(0,5e-324) 或在浏览器的javascript 控制台中输入Math.random()*5e-324。

即使没有上溢/下溢/反规范，也很难可靠地推断浮点操作。经过一番挖掘，我可以找到一个反例：

>>> a=1.0
>>> b=2**-54
>>> rand=a-2*b
>>> a
1.0
>>> b
5.551115123125783e-17
>>> rand
0.9999999999999999
>>> (a-b)*rand+b
1.0

更容易解释为什么 a=2⁵³ 和 b=0.5 会发生这种情况：2⁵³-1 是下一个可表示的数字。默认舍入模式（“舍入到最接近的偶数”）向上舍入 2⁵³-0.5（因为 2⁵³ 是“偶数” [LSB = 0] 和 2 ⁵³-1 是“奇数”[LSB = 1])，所以减去 b 得到 2⁵³，乘以得到 2⁵³-1 ，并添加b 再次获得2⁵³。

---

回答你的第二个问题：因为底层 PRNG 几乎总是在区间 [0,2ⁿ-1] 内生成一个随机数，即它生成随机位。选择一个合适的 n（浮点表示中的精度位）并除以 2ⁿ 并获得可预测的分布非常容易。请注意，[0,1) 中有一些数字，您将永远使用此方法生成（IEEE doubles 中的 (0,2^-53) 中的任何数字）。 p>

这也意味着你可以做到a[Math.floor(Math.random()*a.length)]而不用担心溢出（作业：在IEEE二进制浮点中，证明b < 1意味着a*b < a代表正整数a）。

另一个好处是您可以将每个随机输出 x 视为表示一个区间 [x,x+2^-53)（不太好的事情是平均值返回略小于 0.5)。如果返回 [0,1]，返回端点的概率是否与其他所有端点相同，或者它们是否应该只有一半的概率，因为它们只代表一半的区间？

为了回答在 [0,1] 中返回一个数字的简单问题，下面的方法有效地生成一个整数 [0,2ⁿ]（通过在 [0,2n+1-1]，如果太大就扔掉）除以2ⁿ：

function randominclusive() 
  // Generate a random "top bit". Is it set?
  while (Math.random() >= 0.5) 
    // Generate the rest of the random bits. Are they zero?
    // If so, then we've generated 2^n, and dividing by 2^n gives us 1.
    if (Math.random() == 0)  return 1.0; 
    // If not, generate a new random number.
  
  // If the top bits are not set, just divide by 2^n.
  return Math.random();

cmets 暗示基数为 2，但我认为假设如下：

0 和 1 应该以相等的概率返回（即 Math.random() 不使用接近 0 的浮点数的更近间距）。 Math.random() >= 0.5，概率为 1/2（对于偶数基数应该为真） 底层 PRNG 足够好，我们可以做到这一点。

请注意，随机数总是成对生成：while (a) 中的一个总是跟在 if 中的一个或末尾的一个 (b) 之后。通过考虑返回 0 或 0.5 的 PRNG，很容易验证它是否合理：

a=0   b=0  ：返回0 a=0   b=0.5: 返回 0.5 a=0.5 b=0  ：返回1 a=0.5 b=0.5: 循环

问题：

假设可能不正确。特别是，常见的 PRNG 是采用 48 位 LCG 的前 32 位（Firefox 和 Java 这样做）。要生成双精度数，您需要从两个连续输出中取 53 位并除以 2⁵³，但有些输出是不可能的（您无法生成 2⁵³ 48 位输出国家！）。我怀疑其中一些永远不会返回 0（假设是单线程访问），但我现在不想检查 Java 的实现。 Math.random() 对于每个 潜在 输出都是两次，因为需要获得额外的位，但这对 PRNG 施加了更多限制（需要我们推理四个连续输出上述 LCG）。 Math.random() 每个输出平均被调用 四次 次。有点慢。 它会确定性地丢弃结果（假设单线程访问），因此几乎可以保证减少输出空间。

【参考方案8】：

我对这个问题的解决方案一直是使用以下代替你的上限。

Math.nextAfter(upperBound,upperBound+1)

或

upperBound + Double.MIN_VALUE

所以你的代码应该是这样的：

double myRandomNum = Math.random() * Math.nextAfter(upperBound,upperBound+1) + lowerBound;

或

double myRandomNum = Math.random() * (upperBound + Double.MIN_VALUE) + lowerBound;

这只是将您的上限增加最小的双倍 (Double.MIN_VALUE)，以便您的上限将作为随机计算的可能性包含在内。

这是一个很好的解决方法，因为它不会使概率偏向任何一个数字。

这不起作用的唯一情况是您的上限等于Double.MAX_VALUE

【讨论】：

..以及将Double.MAX_VALUE 作为上限处理的正确代码是什么？

【参考方案9】：

private static double random(double min, double max) 
    final double r = Math.random();
    return (r >= 0.5d ? 1.5d - r : r) * (max - min) + min;

【讨论】：

也许您可以解释一下您的代码中发生了什么，以帮助其他人。

它从结果范围中排除了0.5 值，但包括0 和1。我喜欢这个答案 - ***.com/a/9724775/658606 但该代码两次调用 Math.random。

【参考方案10】：

我的经验相当少，所以我也在寻找解决方案。

这是我的粗略想法：

随机数生成器生成 [0,1) 中的数字，而不是 [0,1]，

因为[0,1)是一个单位长度，后面可以跟[1,2)等不重叠。

对于随机[x, y]， 你可以这样做：

float randomInclusive(x, y)

    float MIN = smallest_value_above_zero;
    float result;
    do
        result = random(x, (y + MIN));
     while(result > y);
    return result;

其中 [x, y] 中的所有值都有相同的被选取的可能性，您现在可以到达 y。

【参考方案11】：

Math.round() 将有助于包含绑定值。如果您有0 <= value < 1（1 是独占），那么Math.round(value * 100) / 100 将返回0 <= value <= 1（1 是包括在内）。这里需要注意的是，该值现在只有 2 位小数位。如果您想要 3 位数字，请尝试 Math.round(value * 1000) / 1000 等等。下面的函数多了一个参数，就是小数位的位数——我称之为precision：

function randomInRange(min, max, precision) 
    return Math.round(Math.random() * Math.pow(10, precision)) /
            Math.pow(10, precision) * (max - min) + min;

【参考方案12】：

这个怎么样？

function randomInRange(min, max)
    var n = Math.random() * (max - min + 0.1) + min;
    return n > max ? randomInRange(min, max) : n;

如果你在这方面遇到堆栈溢出，我会给你买礼物。

-- 编辑：别管现在。我很疯狂：

randomInRange(0, 0.0000000000000000001)

堆栈溢出。

【参考方案13】：

在一个范围内生成一个随机浮点数并非易事。例如，将一个随机整数乘以或除以一个常数，或将一个统一的浮点数缩放到所需范围的常见做法，其缺点是并非浮点格式可以在该范围内表示的所有数字都可以以这种方式覆盖，并且可能存在细微的偏差问题。这些问题在 F. Goualard 的“Generating Random Floating-Point Numbers by Dividing Integers: a Case Study”中有详细讨论。

为了说明问题的重要性，以下伪代码在闭区间 [lo, hi] 中生成一个随机浮点数，其中数字的形式为 FPSign * FPSignificand * FPRADIX^FPExponent。下面的伪代码是从我在floating-point number generation 上的部分复制的。请注意，它适用于浮点数的任何精度和任何基数（包括二进制和十进制）。

METHOD RNDRANGE(lo, hi)
  losgn = FPSign(lo)
  hisgn = FPSign(hi)
  loexp = FPExponent(lo)
  hiexp = FPExponent(hi)
  losig = FPSignificand(lo)
  hisig = FPSignificand(hi)
  if lo > hi: return error
  if losgn == 1 and hisgn == -1: return error
  if losgn == -1 and hisgn == 1
    // Straddles negative and positive ranges
    // NOTE: Changes negative zero to positive
    mabs = max(abs(lo),abs(hi))
    while true
       ret=RNDRANGE(0, mabs)
       neg=RNDINT(1)
       if neg==0: ret=-ret
       if ret>=lo and ret<=hi: return ret
    end
  end
  if lo == hi: return lo
  if losgn == -1
    // Negative range
    return -RNDRANGE(abs(lo), abs(hi))
  end
  // Positive range
  expdiff=hiexp-loexp
  if loexp==hiexp
    // Exponents are the same
    // NOTE: Automatically handles
    // subnormals
    s=RNDINTRANGE(losig, hisig)
    return s*1.0*pow(FPRADIX, loexp)
  end
  while true
    ex=hiexp
    while ex>MINEXP
      v=RNDINTEXC(FPRADIX)
      if v==0: ex=ex-1
      else: break
    end
    s=0
    if ex==MINEXP
      // Has FPPRECISION or fewer digits
      // and so can be normal or subnormal
      s=RNDINTEXC(pow(FPRADIX,FPPRECISION))
    else if FPRADIX != 2
      // Has FPPRECISION digits
      s=RNDINTEXCRANGE(
        pow(FPRADIX,FPPRECISION-1),
        pow(FPRADIX,FPPRECISION))
    else
      // Has FPPRECISION digits (bits), the highest
      // of which is always 1 because it's the
      // only nonzero bit
      sm=pow(FPRADIX,FPPRECISION-1)
      s=RNDINTEXC(sm)+sm
    end
    ret=s*1.0*pow(FPRADIX, ex)
    if ret>=lo and ret<=hi: return ret
  end
END METHOD