生成一个随机的二进制数，其中1位的比例由变量控制。

Question

生成一个随机的二进制数，其中1位的比例由变量控制。

9

我需要一个能够生成随机整数的函数（暂时假设为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值？

- finnw

你好，出于好奇，你测试过被接受的答案是否有效吗？需要多少次迭代才能接近给定的概率？性能仍然足够吗？ - Wim Coenen

是的，我将迭代次数限制在16次（1次迭代等于1位精度），在我的（相当老旧的）个人电脑上每秒生成大约1700万位。我使用了无理概率（例如0.1*PI）进行测量，但对于整数来说速度要快得多，例如p=0.75时每秒生成120兆位。 - finnw

1

有什么特别的原因要避免位移和乘法运算吗？ - Svante

@Svante，不是的。我的第一次尝试太慢了，我以为是因为位移/乘法，但实际上是因为我生成了太多的随机数（通过为每个输出位单独掷骰子）。 - finnw

7个回答

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 );
}

我希望你能理解我的意思。也许可以用long和位运算来替换byte[]，以使其更快。

- Ondra Žižka

我希望我理解了这个问题 - 你需要保证位数为1的数量，即某些位的排列吗？ - Ondra Žižka

有趣，但不符合位必须彼此独立的要求。该方法将始终在结果中返回确切的（64*P）个1。因此，例如当P=1/16时，恰好设置4个位，因此如果设置了位0-3，则所有其他位必须清除。 - finnw

然而，如果将这个方法与一个二项分布生成器结合起来选择“bitCount”，那么它就可以得出正确的结果。 - finnw

我现在将其实现为（a）用于bitCount的二项式分布生成器和（b) 比特位排列的查找表的组合。 - finnw

2

以下是我的最终解决方案：

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

这部分内容受Ondra Žižka答案的启发。

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

当然，计算二项式概率同样耗费时间，因此它们放在另一个查找表中。

虽然代码很多，但在性能方面它正在产生回报。

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

- finnw

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

1

如果您想应用一些分布，其中以概率P获得1，并且以概率1-P获得0，在任何特定位上，那么您最好的选择就是简单地独立生成每个二进制位，并具有成为1的P概率（我知道这听起来像是递归定义）。

以下是解决方案，我将在下面详细介绍：

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个1。

- Mark Elliot

我的原始实现看起来与这个非常相似，但速度非常慢。这就是为什么我决定尝试位操作来并行处理单词中的位。如果P是像0.25、0.5或0.75这样的圆数，它会给出巨大的性能提升。但我不确定对于其他P值是否也适用。 - finnw

1

这是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;
}

- finnw

这看起来很不错。我方法的唯一缺点是你不能指定一个容差在概率上限制生成随机数的数量。但是实现尾数中位数的修剪会比使用tol获得更好的结果，以确保减少运行时间更容易。 - Michael Anderson

实际上，对于非常小的p，内部循环存在问题。我认为-exponent可能比shiftReg的长度更大，在这些情况下，您会做比必要更多的工作。我还担心它可能对非规范化数字做错事情。也许您可以执行long binrep = p * 0xFFFFFFFF并循环遍历binrep中的位，而不是玩弄浮点数的位？然后，您将具有固定长度的循环和近似中已知的最大误差。 - Michael Anderson

可以通过首先调用 iterationCount() 来应用容差，如果答案太高则使用另一种方法。单精度浮点数的最小非零值为 2 ** -149，因此最坏情况下的迭代次数将为 149。异常浮点数在第 9 行进行处理（mantissa |= 0x800000 恢复正常数的隐含 1 位，但跳过了异常浮点数）。如果 -exponent 很大，则不幸的是您需要循环那么多次才能降低概率，即使对于大多数迭代 shiftReg 都为零。 - finnw

但它无法正确处理零。我已经为P = 0添加了一个特殊情况。 - finnw

1

使用一个随机生成器，生成一个在0到1之间的均匀浮点数r。如果r>p，则将位设置为0，否则将其设置为1。

- President James K. Polk

这就是我试图避免的。对于我的应用程序来说太慢了。 - finnw

啊，我明白你的意思了。我会再考虑一下的。 - President James K. Polk

通过将位的批次计算在一起，你不会获得数量级的性能提升。 - Tom Hawtin - tackline

1

@Tom，实际上，在p=0.25的情况下，我确实获得了一个数量级的改进，因为我为每个输出单词生成了4个随机单词（假设RNG输出宽度为32位）。如果我一次只生成一个比特，那么我必须生成64个随机单词，将成本增加了16倍。 - finnw

0

假设位数组的大小为L。如果L=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);
        }
    }
}

EDIT: This code does generate completely random bits. I will try to explain the algorithm better.

Each bit string has a certain probability of occurring. Suppose a string has a probability of occurrence p; we want to choose that string if our random number falls is some interval of length p. The starting point of the interval must be fixed, but its value will not make much difference. Suppose we have chosen upto k bits correctly. Then, for the next bit, we divide the interval corresponding to this k-length bit-string into two parts of sizes in the ratio P:1-P (here P is the probability of getting a 1). We say that the next bit will be 1 if the random number is in the first part, 0 if it is in the second part. This ensure that the probabilities of strings of length k+1 also remain correct.

Java code:

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;
}

- Kartik Kale

@finnw，我看到你没有正确理解算法。我已经改进了解释，请仔细阅读。至于错误的语言，我认为即使对于不熟悉JavaScript的人来说，只要将其视为伪代码，代码也很容易理解。不过，我还是会添加Java代码。 - Kartik Kale

如果double具有无限精度，这将起作用 :-) 很抱歉，但我确实尝试运行了您的代码，并且如我所预期的那样，低阶输出位之间存在很多相关性。 - finnw

是的，只有当您想要与 double 大小相同数量的随机位时，这将起作用。（更准确地说，是随机数生成器的熵，称为 E）。因此，您可以每 E 位使用一个随机数（就像您在解决方案中所做的那样，生成位块）。但是，任何使用固定数量的随机数的解决方案都存在问题。毕竟，您不能从一次抛硬币中挤出多个随机位。除此之外，我不明白为什么您“期望”此解决方案在位之间生成相关性。 - Kartik Kale

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Michael Anderson · Accepted Answer

首先，让我们看一下您已经在代码中使用的一些不太好看的数学。

定义x和y是位，其1的概率分别为X = p（x=1），Y = p（y=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;
}