生成具有可变比例“1”位的随机二进制数答案

【问题标题】：Generate a random binary number with a variable proportion of '1' bits生成具有可变比例“1”位的随机二进制数
【发布时间】：2011-01-05 18:19:19
【问题描述】：

我需要一个函数来生成随机整数。（暂时假设 Java 类型为 long，但稍后将扩展为 BigInteger 或 BitSet。）

棘手的部分是有一个参数 P 指定结果中任何位的（独立）概率为 1。

如果 P = 0.5，那么我们可以使用标准随机数生成器。 P 的其他一些值也很容易实现。这是一个不完整的例子：

Random random = new Random();

// ...

long nextLong(float p) {
    if      (p == 0.0f)   return 0L;
    else if (p == 1.0f)   return -1L;
    else if (p == 0.5f)   return random.nextLong();
    else if (p == 0.25f)  return nextLong(0.5f) & nextLong(0.5f);
    else if (p == 0.75f)  return nextLong(0.5f) | nextLong(0.5f);
    else if (p == 0.375f) return nextLong(0.5f) & nextLong(0.75f); // etc
    else {
      // What goes here??
      String message = String.format("P=%f not implemented yet!", p);
      throw new IllegalArgumentException(message);
    }
}

对于介于 0.0 和 1.0 之间的任何 P 值，有没有办法对此进行概括？

【问题讨论】：

出于好奇，您是否测试过接受的答案是否有效？需要多少次迭代才能可接受地接近给定的概率？性能还不够吗？
是的，我将其限制为 16 次迭代（1 次迭代 == 1 位精度），它在我的（相当旧的）PC 上每秒生成大约 1700 万位。我用非理性概率（例如 0.1*PI）测量了这个，但是对于整数来说它要快得多，例如p=0.75 时为 120Mbits/秒。
有什么特别的理由可以避免移位和乘法吗？
@Svante，没有。我的第一次尝试太慢了，我认为这是因为位移/乘法，但实际上只是我生成了太多随机数（通过为每个输出位单独滚动。）

标签： java optimization random bit-manipulation

【解决方案1】：

首先是您已经在代码中使用的有点难看的数学。

定义 x 和 y 分别是 X = p(x=1), Y = p(y=1) 中概率为 1 的位。然后我们就有了

 p( x & y = 1) = X Y
 p( x | y = 1) = 1 - (1-X) (1-Y)
 p( x ^ y = 1) = X (1 - Y) + Y (1 - X)

现在如果我们让 Y = 1/2 我们得到

P( x & y ) = X/2
P( x | y ) = (X+1)/2

现在将 RHS 设置为我们想要的概率，我们有两种情况可以解决 X

X = 2 p        // if we use &
X = 2 p - 1    // if we use |

接下来我们假设我们可以再次使用它来根据另一个变量 Z 获得 X... 然后我们继续迭代，直到我们完成“足够”。

这有点不清楚，但考虑 p = 0.375

0.375 * 2 = 0.75  < 1.0 so our first operation is &
0.75 * 2 = 1.5 > 1.0 so our second operation is |
0.5 is something we know so we stop.

因此我们可以通过 X1 & (X2 | X3) 得到一个 p=0.375 的变量

问题是对于大多数变量，这不会终止。例如

0.333 *2 = 0.666 < 1.0 so our first operation is &
0.666 *2 = 1.333 > 1.0 so our second operation is |
0.333 *2 = 0.666 < 1.0 so our third operation is &
etc...

所以p=0.333可以由

生成

X1 & ( X2 | (X3 & (X4 | ( ... ) ) ) )

现在我怀疑在系列中取足够多的术语会给你足够的准确性，这可以写成递归函数。然而，也可能有更好的方法......我认为操作的顺序与 p 的二进制表示有关，我只是不确定如何......并且没有时间更深入地考虑它。

无论如何，这里有一些未经测试的 C++ 代码可以做到这一点。您应该能够轻松地对其进行 java 化。

uint bitsWithProbability( float p )
{
   return bitsWithProbabilityHelper( p, 0.001, 0, 10 );
}

uint bitsWithProbabilityHelper( float p, float tol, int cur_depth, int max_depth )
{
   uint X = randbits();
   if( cur_depth >= max_depth) return X;
   if( p<0.5-tol)
   {
     return X & bitsWithProbabilityHelper( 2*p, 0.001, cur_depth+1, max_depth );
   }
   if(p>0.5+tol)
   {
     return X | bitsWithProbabilityHelper( 2*p-1, 0.001, cur_depth+1, max_depth );
   }
   return X;
}

【讨论】：

您可以在每个步骤调整 tol 并删除 max_depth。再次只需要比我目前拥有的更多的备用大脑。
+1 我曾经做过同样的事情，正要测试它是否收敛
是的，我认为cur_depth 是多余的。将tol 乘以2.001 等常数因子会产生类似的效果。
@finw 确实。事实上，在 N 项之后截断序列的概率误差最多为 pow(2,-N)。因此，每次迭代将 tol 缩放 2 是明智/正确的。
我想出了一种使用 P 的二进制表示进行迭代的方法。请参阅我的第二个答案。

【解决方案2】：

在数字中按比例分配位数。伪代码：

long generateNumber( double probability ){
  int bitCount = 64 * probability;
  byte[] data = new byte[64]; // 0-filled

  long indexes = getRandomLong();

  for 0 to bitCount-1 {
    do { 
      // distribute this bit to some postition with 0.
      int index = indexes & 64;
      indexes >> 6;
      if( indexes == 0 ) indexes = getRandomLong();
    } while ( data[index] == 0 );
    data[index] = 1;
  }

  return bytesToLong( data );
}

我希望你明白我的意思。也许byte[] 可以替换为long 和位操作以使其更快。

【讨论】：

我希望我理解了这个问题 - 您是否需要保证位数为 1，即某些位排列？
有趣，但不满足位必须相互独立的要求。此方法将始终在结果中准确返回 (64*P) 个。因此，例如 P=1/16 恰好设置了 4 位，因此如果设置了 0-3 位，则所有其他位都必须清除。
但是，如果将其与二项分布生成器结合使用以选择 bitCount，那么它可以给出正确的结果。
我现在已经将其实现为 (a) bitCount 的二项分布生成器和 (b) 位排列查找表的组合。

【解决方案3】：

这就是我最终解决的方法。

按照二项分布生成一个介于 0..16 之间的整数 N。这给出了 16 位部分结果中“1”的位数。
在包含 16 位整数的查找表中随机生成一个索引，其中包含所需数量的“1”位。
重复 4 次，得到四个 16 位整数。
将这四个 16 位整数拼接在一起，得到一个 64 位整数。

这部分灵感来自 Ondra Žižka 的回答。

好处是它将对Random.nextLong() 的调用次数减少到每 64 位输出调用 8 次。相比之下，滚动每个单独的位将需要 64 次调用。根据 P 的值，按位 AND/OR 使用 2 到 32 次调用

当然，计算二项式概率的成本同样高，所以它们会放在另一个查找表中。

这是很多代码，但它在性能方面得到了回报。

更新 - 将其与按位 AND/OR 解决方案合并。如果它猜测它会更有效（就对Random.next()的调用而言），它现在会使用该方法。

【讨论】：

我也在考虑这个问题，并且独立想出了类似你的按位与/或方法。我可以检查我对您的最终解决方案的理解吗？您有两个 LUT：1) 包含 n=16 和 P=0.5 的二项式分布，以及 2) 包含按设置的位数分组的所有 2^16 16 位整数。您生成一个随机索引以从 LUT 1 中选择多个位，然后生成另一个随机索引以从 LUT 2 中选择具有那么多位的 16 位值。这样正确吗？

【解决方案4】：

使用随机生成器生成介于 0 和 1 之间的统一浮点数 r。如果 r>p 则将该位设置为 0，否则将其设置为 1

【讨论】：

这就是我要避免的。我的应用程序太慢了。
啊，我明白你在说什么。我会考虑更多。
通过一起计算成批的位，您不会获得数量级的性能提升。
@Tom，实际上我确实得到了数量级的改进，至少在 p=0.25 的情况下，因为我为每个输出词生成 4 个随机词（假设 RNG 输出是32 位宽。）如果我一次生成一位，我必须生成 64 个随机字，成本增加了 16 倍。

【解决方案5】：

如果您要应用某种分布，其中概率 P 得到 1，概率 1-P 在任何特定位得到 0，您最好的选择就是独立生成每个位，概率为 P 1（我知道这听起来像是一个递归定义）。

这是一个解决方案，我将在下面逐步介绍：

public class MyRandomBitGenerator
{

    Random pgen = new Random();

    // assumed p is well conditioned (0 < p < 1)
    public boolean nextBitIsOne(double p){
        return pgen.nextDouble() < p ? true : false;
    }

    // assumed p is well conditioned (0 < p < 1)
    public long nextLong(double p){
        long nxt = 0;
        for(int i = 0; i < 64; i++){
           if(nextBitIsOne(p)){
               nxt += 1 << i;
           }
        }
        return nxt;
    }

}

基本上，我们首先确定如何以概率 P 生成值 1：pgen.nextDouble() 以均匀概率生成介于 0 和 1 之间的数字，通过询问它是否小于 p，我们正在对该分布进行采样，使得我们希望看到p 1s，因为我们无限地调用这个函数。

【讨论】：

我的原始实现看起来与此非常相似，而且非常慢。这就是为什么我决定尝试位操作来并行化一个单词中的位。如果 P 是 0.25、0.5 或 0.75 之类的整数，它会带来巨大的性能提升。我不确定这是否适用于其他 P 值。

【解决方案6】：

这是Michael Anderson's answer的另一个变体

为了避免递归，我们从右到左迭代处理 P 的位，而不是从左到右递归处理。这在浮点表示中会很棘手，因此我们从二进制表示中提取指数/尾数字段。

class BitsWithProbabilityHelper {
    public BitsWithProbabilityHelper(float prob, Random rnd) {
        if (Float.isNaN(prob)) throw new IllegalArgumentException();

        this.rnd = rnd;

        if (prob <= 0f) {
            zero = true;
            return;
        }

        // Decode IEEE float
        int probBits = Float.floatToIntBits(prob);
        mantissa = probBits & 0x7FFFFF;
        exponent = probBits >>> 23;

        // Restore the implicit leading 1 (except for denormals)
        if (exponent > 0) mantissa |= 0x800000;
        exponent -= 150;

        // Force mantissa to be odd
        int ntz = Integer.numberOfTrailingZeros(mantissa);
        mantissa >>= ntz;
        exponent += ntz;
    }

    /** Determine how many random words we need from the system RNG to
     *  generate one output word with probability P.
     **/
    public int iterationCount() {
        return - exponent;
    }

    /** Generate a random number with the desired probability */
    public long nextLong() {
        if (zero) return 0L;

        long acc = -1L;
        int shiftReg = mantissa - 1;
        for (int bit = exponent; bit < 0; ++ bit) {
            if ((shiftReg & 1) == 0) {
                acc &= rnd.nextLong();
            } else {
                acc |= rnd.nextLong();
            }
            shiftReg >>= 1;
        }
        return acc;
    }

    /** Value of <code>prob</code>, represented as m * 2**e where m is always odd. */
    private int exponent;  
    private int mantissa;

    /** Random data source */
    private final Random rnd;

    /** Zero flag (special case) */
    private boolean zero;
}

【讨论】：

这看起来不错。我的方法的唯一缺点是您无法指定概率的容差来限制生成的随机数的数量。但是，很容易实现尾数中的位剪辑，这将比使用 tol 提供更好的结果，从而保证减少运行时间。
实际上，对于非常小的 p，内部循环有问题。我认为 -exponent 可以大于 shiftReg 的长度，在这种情况下，您所做的工作超出了必要的范围。我还担心它可能对去正常化的数字做错了。也许你可以做一个 long binrep = p * 0xFFFFFFFF 并循环遍历 binrep 中的位而不是玩浮点数？然后你会有一个固定长度的循环和近似值中已知的最大误差。
您可以先调用iterationCount() 来应用容差，如果答案太高，请使用其他方法。单精度浮点数的最小非零值为 2**-149，因此最坏情况下的迭代计数为 149。非规范化数字在第 9 行处理（mantissa |= 0x800000 恢复正常数字的隐式 1 位，对于非规范化的数字，这将被跳过。）如果 -exponent 很大，那么不幸的是，您确实需要循环多次以降低概率，即使对于大多数迭代来说 shiftReg 将为零。
它没有正确处理零。我为 P=0 添加了一个特殊情况。

【解决方案7】：

假设位数组的大小为L。如果L=1，则第1位为1的概率为P，为0的概率为1-P。对于 L=2，得到 00 的概率是 (1-P)²，01 或 10 分别是 P(1-P)，11 是 P² .扩展这个逻辑，我们可以首先通过将随机数与 P 进行比较来确定第一位，然后缩放随机数，这样我们就可以再次得到 0 到 1 之间的任何值。示例 javascript 代码：

function getRandomBitArray(maxBits,probabilityOf1) {
    var randomSeed = Math.random();
    bitArray = new Array();
    for(var currentBit=0;currentBit<maxBits;currentBit++){
        if(randomSeed<probabilityOf1){
            //fill 0 at current bit
            bitArray.push(0);
            //scale the sample space of the random no from [0,1)
            //to [0.probabilityOf1)
            randomSeed=randomSeed/probabilityOf1;
        }
        else{
            //fill 1 at current bit
            bitArray.push(1);
            //scale the sample space to [probabilityOf1,1)
            randomSeed = (randomSeed-probabilityOf1)/(1-probabilityOf1);
        }
    }
}

编辑： 此代码确实会生成完全随机的位。我会尝试更好地解释算法。

每个位串都有一定的发生概率。假设一个字符串有一个出现的概率p；如果我们的随机数落在某个长度为 p 的区间，我们想选择那个字符串。区间的起点必须是固定的，但它的值不会有太大的不同。假设我们正确选择了 k 位。然后，对于下一个比特，我们将这个 k 长度比特串对应的区间分成大小比例为 P:1-P 的两部分（这里 P 是得到 1) 的概率。如果随机数在第一部分，我们说下一位为 1，如果在第二部分，则为 0。这确保了长度为 k+1 的字符串的概率也保持正确。

Java 代码：

public ArrayList<Boolean> getRandomBitArray(int maxBits, double probabilityOf1) {
    double randomSeed = Math.random();
    ArrayList<Boolean> bitArray = new ArrayList<Boolean>();
    for(int currentBit=0;currentBit<maxBits;currentBit++){
        if(randomSeed<probabilityOf1){
            //fill 0 at current bit
            bitArray.add(false);
            //scale the sample space of the random no from [0,1)
            //to [0.probabilityOf1)
            randomSeed=randomSeed/probabilityOf1;
        }
        else{
            //fill 1 at current bit
            bitArray.add(true);
            //scale the sample space to [probabilityOf1,1)
            randomSeed = (randomSeed-probabilityOf1)/(1-probabilityOf1);
        }
    }
    return  bitArray;
}

【讨论】：

@finnw 我看到你没有正确理解算法。我改进了解释。请仔细阅读。关于错误的语言，我认为如果将代码视为伪代码，即使对于不熟悉 javascript 的人来说，代码也很容易理解。不过，我正在添加 Java 代码。
如果只有double 具有无限精度，这将起作用:-) 抱歉，我确实尝试过运行您的代码，并且正如我所料，低位输出位之间存在很多相关性。
是的，只有当您只需要与 double 大小一样多的随机位时，这才有效。（更准确地说，是随机数生成器的熵，比如 E）。因此，您可以为每个 E 位使用一个随机数（就像您在解决方案中所做的那样，生成位块）。但这对于使用固定数量的随机数的任何解决方案都是一个问题。毕竟，你不能以某种方式从抛硬币中挤出一个以上的随机比特。除此之外，我不明白您为什么“期望”此解决方案会在位之间产生相关性。