快速生成伪随机位的方法，每个位的给定概率为 0 或 1
Fast way to generate pseudo-random bits with a given probability of 0 or 1 for each bit
通常，随机数生成器 returns 比特流，在每个位置观察到 0 或 1 的概率是相等的（即 50%）。我们称其为无偏 PRNG。
我需要生成一串具有以下内容的伪随机位 属性：在每个位置看到 1 的概率是 p（即看到 0 的概率是 1-p）。参数p是0到1之间的实数；在我的问题中，碰巧它的分辨率为 0.5%，即它可以取值 0%、0.5%、1%、1.5%、...、99.5%、100%。
请注意，p 是概率而不是精确分数。 n位流中实际设置为1的位数必须服从二项分布B(n, p).
有一种天真的方法可以使用无偏PRNG来生成每个位的值（伪代码）：
generate_biased_stream(n, p):
result = []
for i in 1 to n:
if random_uniform(0, 1) < p:
result.append(1)
else:
result.append(0)
return result
这样的实现比生成无偏流的实现要慢得多，因为它为每个位调用一次随机数生成器函数；而无偏流生成器每个字大小调用一次（例如，它可以通过一次调用生成 32 或 64 个随机位）。
我想要一个更快的实现，即使它稍微牺牲了随机性。想到的一个想法是预先计算查找 table：对于 p 的 200 个可能值中的每一个，使用较慢的算法计算 C 8 位值并将它们保存在 table 中。然后快速算法会随机选择其中一个来生成 8 个倾斜位。
用于查看需要多少内存的粗略计算：
C 至少应为 256（可能的 8 位值的数量），可能更多以避免采样效应；比方说 1024。也许这个数字应该根据 p 而变化，但让我们保持简单并假设平均值是 1024。
由于 p => 有 200 个值 => 总内存使用量为 200 KB。这还不错，可能适合 L2 缓存 (256 KB)。我还需要评估一下，看看是否有抽样效应引入偏差，这种情况下C就得调高了。
此解决方案的不足之处在于它一次只能生成 8 位，即使需要大量工作，而无偏 PRNG 只需几条算术指令即可一次生成 64 位。
我想知道是否有更快的方法，基于位操作而不是查找 tables。例如直接修改随机数生成代码，为每一位引入偏差。这将实现与无偏 PRNG 相同的性能。
编辑 3 月 5 日
谢谢大家的建议，我得到了很多有趣的想法和建议。以下是最重要的：
- 更改问题要求，使 p 的分辨率为 1/256 而不是 1/200。这允许更有效地使用位，并且还提供了更多的优化机会。我想我可以做出这个改变。
- 使用算术编码有效地消耗来自无偏生成器的位。通过上述分辨率更改，这变得容易多了。
- 一些人认为 PRNG 非常快，因此使用算术编码实际上可能会由于引入的开销而使代码变慢。相反，我应该始终使用最坏情况下的位数并优化该代码。请参阅下面的基准。
- @rici 建议使用 SIMD。这是一个不错的想法，只有当我们总是消耗固定数量的比特时它才有效。
基准（没有算术解码）
注意：正如你们许多人所建议的，我将分辨率从 1/200 更改为 1/256。
我写了 several implementations 简单地采用 8 个随机无偏位并生成 1 个有偏位的朴素方法：
- 没有 SIMD
- 使用 Agner Fog 的向量类库的 SIMD，
- 使用内部函数的 SIMD
我使用两个无偏伪随机数生成器：
- xorshift128plus
- Ranvec1 (Mersenne Twister-like) 来自 Agner Fog 的图书馆。
我还测量了无偏 PRNG 的速度以进行比较。以下是结果：
RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 16.081 16.125 16.093 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 0.778 0.783 0.812 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.176 2.184 2.145 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.129 2.151 2.183 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
与标量方法相比，SIMD 将性能提高了 3 倍。正如预期的那样，它比无偏生成器慢 8 倍。
最快的偏置生成器达到 2.1 Gb/s。
RNG: xorshift128plus
Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.300 21.486 21.483 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical : 104,857,600
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 22.660 22.661 24.662 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 1.065 1.102 1.078 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 4.972 4.971 4.970 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 4.955 4.971 4.971 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical : 104,857,600
对于 xorshift，与标量方法相比，SIMD 将性能提高了 5 倍。它比无偏生成器慢 4 倍。请注意，这是 xorshift 的标量实现。
最快的偏置生成器达到 4.9 Gb/s。
RNG: xorshift128plus_avx2
Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.754 21.494 21.878 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical : 104,857,600
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 54.126 54.071 54.145 [Gb/s]
Number of ones: 536,874,540 536,880,718 536,891,316
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 1.093 1.103 1.063 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 19.567 19.578 19.555 [Gb/s]
Number of ones: 104,836,115 104,846,215 104,835,129
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 19.551 19.589 19.557 [Gb/s]
Number of ones: 104,831,396 104,837,429 104,851,100
Theoretical : 104,857,600
此实现将 AVX2 用于 运行 4 个并行的无偏异或移位生成器。
最快的偏置生成器达到 19.5 Gb/s。
算术解码基准
简单的测试表明算术解码代码是瓶颈，而不是PRNG。所以我只对最昂贵的 PRNG 进行基准测试。
RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)
Method: Arithmetic decoding (floating point)
Gbps/s: 0.068 0.068 0.069 [Gb/s]
Number of ones: 10,235,580 10,235,580 10,235,580
Theoretical : 10,240,000
Method: Arithmetic decoding (fixed point)
Gbps/s: 0.263 0.263 0.263 [Gb/s]
Number of ones: 10,239,367 10,239,367 10,239,367
Theoretical : 10,240,000
Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 12.687 12.686 12.684 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical : 104,857,600
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 14.536 14.536 14.536 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 0.754 0.754 0.754 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.094 2.095 2.094 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.094 2.094 2.095 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
简单的不动点方法达到 0.25 Gb/s，而朴素的标量方法快 3 倍，朴素的 SIMD 方法快 8 倍。可能有一些方法可以进一步优化 and/or 并行化算术解码方法，但由于其复杂性，我决定就此打住并选择简单的 SIMD 实现。
谢谢大家的帮助。
如果p接近于0，则可以计算出第n位是第一个为1的位的概率；然后你计算一个介于 0 和 1 之间的随机数，并相应地选择 n。例如，如果 p = 0.005 (0.5%)，并且随机数是 0.638128，您可能会计算（我在这里猜测）n = 321，因此您填充 321 个 0 位和一个位集。
如果p接近1，用1-p代替p，设置1位加一个0位。
如果p不接近1或0，对所有256个8位序列做一个table，计算它们的累积概率，然后得到一个随机数，在数组中进行二分查找累积概率，可以设置8位。
假设 1 出现的概率是 6.25% (1/16)。 4 位数字有 16 种可能的位模式：
0000,0001, ..., 1110,1111.
现在，只需像以前一样生成一个随机数，并将半字节边界处的每个 1111 替换为 1，并将其他所有内容变为 0。
根据其他概率进行相应调整。
你可以做的一件事是多次从底层无偏生成器中采样，得到几个32位或64位的字，然后进行按位布尔运算。例如，对于 4 个单词 b1,b2,b3,b4，您可以得到以下分布：
expression | p(bit is 1)
-----------------------+-------------
b1 & b2 & b3 & b4 | 6.25%
b1 & b2 & b3 | 12.50%
b1 & b2 & (b3 | b4) | 18.75%
b1 & b2 | 25.00%
b1 & (b2 | (b3 & b4)) | 31.25%
b1 & (b2 | b3) | 37.50%
b1 & (b2 | b3 | b4)) | 43.75%
b1 | 50.00%
类似的构造可以用于更精细的分辨率。它变得有点乏味并且仍然需要更多的生成器调用，但至少不是每位一个。这类似于 a3f 的答案，但可能更容易实现，而且我怀疑比扫描 0xF nybbles 的单词更快。
请注意，对于您想要的 0.5% 的分辨率，您需要 8 个无偏见的词来对应一个有偏见的词，这将为您提供 (0.5^8) = 0.390625% 的分辨率。
给出精确结果的一种方法是首先为一个 k 位块随机生成符合二项分布的 1 位的数量，然后使用以下方法之一生成具有恰好那么多位的 k 位字方法 here。例如mic006的方法只需要大约log k k位随机数，而我的只需要一个。
它的作用
此实现通过“/dev/urandom”特殊字符文件的接口对随机设备内核模块进行单次调用，以获取表示所有值所需的随机数据数给定的决议。最大可能分辨率为 1/256^2，因此 0.005 可以表示为：
328/256^2,
即：
分辨率：256*256
x: 328
错误 0.000004883。
它是如何做到的
该实现计算位数 bits_per_byte，这是处理给定分辨率所需的均匀分布位数，即表示所有 @resolution 值。然后它会调用随机化设备（如果定义了 URANDOM_DEVICE，则调用“/dev/urandom”，否则它将通过调用“/dev/random”来使用来自设备驱动程序的额外噪声，如果存在，则可能会阻塞没有足够的比特熵）来获得所需数量的均匀分布的字节并填充随机字节数组 rnd_bytes。最后，它从 rnd_bytes 数组的每个 bytes_per_byte 字节中读取每个伯努利样本所需的位数，并将这些位的整数值与 x/resolution 给出的单个伯努利结果的成功概率进行比较。如果值命中，即它落在我们任意选择为 [0, x/resolution) 段的 x/resolution 长度段中，那么我们注意到成功并将 1 插入到结果数组中。
从随机设备读取：
/* if defined use /dev/urandom (will not block),
* if not defined use /dev/random (may block)*/
#define URANDOM_DEVICE 1
/*
* @brief Read @outlen bytes from random device
* to array @out.
*/
int
get_random_samples(char *out, size_t outlen)
{
ssize_t res;
#ifdef URANDOM_DEVICE
int fd = open("/dev/urandom", O_RDONLY);
if (fd == -1) return -1;
res = read(fd, out, outlen);
if (res < 0) {
close(fd);
return -2;
}
#else
size_t read_n;
int fd = open("/dev/random", O_RDONLY);
if (fd == -1) return -1;
read_n = 0;
while (read_n < outlen) {
res = read(fd, out + read_n, outlen - read_n);
if (res < 0) {
close(fd);
return -3;
}
read_n += res;
}
#endif /* URANDOM_DEVICE */
close(fd);
return 0;
}
填写伯努利样本向量：
/*
* @brief Draw vector of Bernoulli samples.
* @details @x and @resolution determines probability
* of success in Bernoulli distribution
* and accuracy of results: p = x/resolution.
* @param resolution: number of segments per sample of output array
* as power of 2: max resolution supported is 2^24=16777216
* @param x: determines used probability, x = [0, resolution - 1]
* @param n: number of samples in result vector
*/
int
get_bernoulli_samples(char *out, uint32_t n, uint32_t resolution, uint32_t x)
{
int res;
size_t i, j;
uint32_t bytes_per_byte, word;
unsigned char *rnd_bytes;
uint32_t uniform_byte;
uint8_t bits_per_byte;
if (out == NULL || n == 0 || resolution == 0 || x > (resolution - 1))
return -1;
bits_per_byte = log_int(resolution);
bytes_per_byte = bits_per_byte / BITS_PER_BYTE +
(bits_per_byte % BITS_PER_BYTE ? 1 : 0);
rnd_bytes = malloc(n * bytes_per_byte);
if (rnd_bytes == NULL)
return -2;
res = get_random_samples(rnd_bytes, n * bytes_per_byte);
if (res < 0)
{
free(rnd_bytes);
return -3;
}
i = 0;
while (i < n)
{
/* get Bernoulli sample */
/* read byte */
j = 0;
word = 0;
while (j < bytes_per_byte)
{
word |= (rnd_bytes[i * bytes_per_byte + j] << (BITS_PER_BYTE * j));
++j;
}
uniform_byte = word & ((1u << bits_per_byte) - 1);
/* decision */
if (uniform_byte < x)
out[i] = 1;
else
out[i] = 0;
++i;
}
free(rnd_bytes);
return 0;
}
用法：
int
main(void)
{
int res;
char c[256];
res = get_bernoulli_samples(c, sizeof(c), 256*256, 328); /* 328/(256^2) = 0.0050 */
if (res < 0) return -1;
return 0;
}
如果您准备基于 256 个可能的值来近似 p，并且您有一个 PRNG 可以生成各个位彼此独立的统一值，那么您可以使用矢量化比较从单个随机数中产生多个偏置位。
只有在 (1) 您担心随机数质量和 (2) 您可能需要大量具有相同偏差的位时才值得这样做。第二个要求似乎是由原始问题暗示的，它批评了一个提出的解决方案，如下： "A deficiency of this solution is that it can generate only 8 bits at once, even that with a lot of work, while an unbiased PRNG can generate 64 at once with just a few arithmetic instructions." 这里，暗示似乎是 useful 来生成一个一次调用中的大块偏置位。
随机数质量是一个困难的课题。即使不是不可能，也很难衡量，因此不同的人会提出不同的衡量标准，强调 and/or 贬低 "randomness" 的不同方面。通常可以权衡随机数生成的速度以获得较低的 "quality"；这是否值得做取决于您的具体应用。
最简单的随机数质量测试涉及单个值的分布和生成器的周期长度。 C 库 rand 和 Posix random 函数的标准实现通常会通过分布测试，但循环长度对于长 运行 应用程序来说是不够的。
不过，这些生成器通常非常快：random 的 glibc 实现只需要几个周期，而经典的线性同余生成器 (LCG) 需要乘法和加法。 （或者，在 glibc 实现的情况下，上面的三个生成 31 位。）如果这足以满足您的质量要求，那么尝试优化就没有意义了，尤其是在偏差概率频繁变化的情况下。
请记住，循环长度应比预期的样本数长很多；理想情况下，它应该大于该数字的平方，因此如果您希望生成千兆字节的随机数，则循环长度为 231 的线性同余生成器 (LCG) 是不合适的数据。即使是声称周期长度约为 235 的 Gnu 三项式非线性加性反馈发生器，也不应用于需要数百万样本的应用程序。
另一个更难测试的质量问题与连续样本的独立性有关。短周期长度在这个指标上完全失败，因为一旦重复开始，生成的随机数就会与历史值精确相关。 Gnu三项式算法虽然周期较长，但由于生成的ith个随机数，ith有明显的相关性=43=]ri,永远是两个值r之一i−3+ri−31 或 ri−3+ri−31+1。这可能会产生 surprising or at least puzzling 后果，尤其是对于伯努利实验。
这是一个使用 Agner Fog 的有用 vector class library, which abstracts away a lot of the annoying details in SSE intrinsics, and also helpfully comes with a fast vectorized random number generator（在 vectorclass.zip 存档中的 special.zip 中找到）的实现，它使我们能够从对 256 位 PRNG 的八次调用中生成 256 位.您可以阅读 Fog 博士关于他为何发现甚至 Mersenne twister 都存在质量问题的解释，以及他提出的解决方案；我真的没有资格发表评论，但它至少在我尝试过的伯努利实验中确实给出了预期的结果。
#include "vectorclass/vectorclass.h"
#include "vectorclass/ranvec1.h"
class BiasedBits {
public:
// Default constructor, seeded with fixed values
BiasedBits() : BiasedBits(1) {}
// Seed with a single seed; other possibilities exist.
BiasedBits(int seed) : rng(3) { rng.init(seed); }
// Generate 256 random bits, each with probability `p/256` of being 1.
Vec8ui random256(unsigned p) {
if (p >= 256) return Vec8ui{ 0xFFFFFFFF };
Vec32c output{ 0 };
Vec32c threshold{ 127 - p };
for (int i = 0; i < 8; ++i) {
output += output;
output -= Vec32c(Vec32c(rng.uniform256()) > threshold);
}
return Vec8ui(output);
}
private:
Ranvec1 rng;
};
在我的测试中，它在 260 毫秒内产生并计数了 268435456 位，或每纳秒一位。测试机是i5，所以没有AVX2； YMMV.
在实际用例中，p 有 201 个可能值，8 位阈值的计算将非常不精确。如果不需要这种不精确性，您可以调整以上内容以使用 16 位阈值，但代价是生成两倍多的随机数。
或者，您可以基于 10 位阈值手动进行矢量化，这将使您非常接近 0.5% 的增量，使用标准的位操作技巧进行矢量化阈值比较，方法是检查借用值向量与重复阈值相减的每第 10 位。与 std::mt19937_64 相结合，每个 64 位随机数平均有 6 位。
假设你可以访问一个随机位生成器，你可以生成一个值来逐位与p进行比较，并在你能证明生成的值小于时立即中止或大于或等于 p.
按照以下步骤以给定的概率在流中创建一个项目 p:
- 以二进制形式
0. 开头
- 附加一个随机位；假设绘制了
1，您将得到 0.1
- 如果结果（二进制表示法）可证明小于
p，则输出1
- 如果结果可证明大于或等于
p，输出一个0
- 否则（如果不能排除两者），继续步骤2。
假设p二进制表示为0.1001101...；如果此过程生成0.0、0.1000、0.10010、...中的任何一个，则该值不能再大于或等于p；如果生成0.11、0.101、0.100111、...中的任何一个，则该值不能小于p。
对我来说，看起来这种方法预期使用了大约两个随机位。对于固定 p，算术编码（如 Mark Dickinson 的回答所示）最多 每个偏置位（平均）消耗一个随机位；修改的成本 p 不清楚。
呃，伪随机数生成器一般都很快。我不确定这是什么语言（Python，也许），但是 "result.append"（几乎肯定包含内存分配）可能比 "random_uniform"（只做一点数学运算）慢.
如果你想优化这段代码的性能：
- 验证这是一个问题。优化需要一些工作，并且会使代码更难维护。除非必要，否则不要这样做。
- 剖析它。 运行 一些测试以确定代码的哪些部分实际上是最慢的。这些是您需要加速的部分。
- 进行更改，并验证它们是否确实更快。编译器非常聪明；通常清晰的代码会编译成更好的代码，比看起来更快的复杂代码。
如果您使用的是编译语言（即使是 JIT 编译语言），每次控制转移（if、while、函数调用等）都会影响性能。消除你能消除的。内存分配也（通常）非常昂贵。
如果您使用的是解释性语言，那么一切皆有可能。最简单的代码很可能是最好的。解释器的开销会让你做什么都相形见绌，所以尽可能减少它的工作。
我只能猜测你的性能问题出在哪里：
- 内存分配。以完整大小预分配数组，稍后填充条目。这确保在添加条目时不需要重新分配内存。
- 分支机构。您可以通过转换结果或类似的东西来避免 "if" 。这在很大程度上取决于编译器。检查程序集（或配置文件）以验证它是否符合您的要求。
- 数值类型。找出您的随机数生成器本机使用的类型，并在该类型中进行算术运算。例如，如果生成器自然 returns 32 位无符号整数，则先将 "p" 缩放到该范围，然后将其用于比较。
顺便说一下，如果您真的想使用尽可能少的随机数，请使用 "arithmetic coding" 来解码您的随机流。不会很快。
您将获得理论上的最佳行为，即真正地最小化使用随机数生成器并能够对任何概率建模p 完全正确， 如果你使用 arithmetic coding.
算术编码是一种数据压缩形式，它将消息表示为数字范围的子区间。它提供了理论上最优的编码，并且可以为每个输入符号使用小数位数。
这个想法是这样的：假设您有一个随机位序列，它们是 1，概率为 p。为方便起见，我将使用 q 作为该位为零的概率。 (q = 1-p)。算术编码分配给数字范围的每个位部分。对于第一位，如果输入为 0，则分配区间 [0, q)，如果输入为 1，则分配区间 [q, 1)。后续位分配当前范围的比例子区间。例如，假设 q = 1/3 输入 1 0 0 将被编码成这样：
Initially [0, 1), range = 1
After 1 [0.333, 1), range = 0.6666
After 0 [0.333, 0.5555), range = 0.2222
After 0 [0.333, 0.407407), range = 0.074074
第一个数字，1，选择范围的前三分之二（1-q）；第二个数字 0 选择 that 的后三分之一，依此类推。
在第一步和第二步之后，区间跨越中点；但是第三步之后完全在中点以下，所以可以输出第一个压缩数字：0。该过程继续，并添加一个特殊的 EOF 符号作为终止符。
这与您的问题有什么关系？压缩后的输出将具有随机的零和等概率的一。因此，要获得概率为 p 的位，只需假设你的 RNG 的输出是上述算术编码的结果，然后 对其应用解码器过程。 也就是说，读取位就好像它们将行间隔细分成越来越小的部分一样。例如，我们从 RNG 中读取 01 后，我们将在 [0.25, 0.5) 范围内。继续读取位，直到足够的输出为 "decoded"。由于您正在模拟解压缩，您将得到比输入更多的随机位。因为算术编码在理论上是最优的，所以无法将 RNG 输出转换为更多有偏差的位在不牺牲随机性的情况下：您获得了真正的最大值。
要注意的是，您无法通过几行代码完成此操作，而且我不知道有什么库可以为您指明（尽管一定有一些您可以使用）。不过，这还是很简单的。 above article provides code for a general-purpose encoder and decoder, in C. It's pretty straightforward, and it supports multiple input symbols with arbitrary probabilities; in your case a far simpler implementation is possible (as 现在显示），因为概率模型是微不足道的。为了扩展使用，需要做更多的工作来产生一个不会对每个位进行大量浮点计算的健壮的实现。
Wikipedia 也有一个有趣的讨论算术编码被认为是基数的变化，这是查看任务的另一种方式。
从信息论的角度来看，有偏差的比特流（p != 0.5）比无偏差的比特流具有更少的信息，所以理论上它应该（平均）少于 比无偏输入的 1 位来产生有偏输出流的一位。例如，p = 0.1的伯努利随机变量的entropy是-0.1 * log2(0.1) - 0.9 * log2(0.9)位，也就是0.469位左右。这表明对于 p = 0.1 的情况，我们应该能够为每个无偏输入位生成两位多一点的输出流。
下面，我给出两种产生偏置位的方法。在需要尽可能少的输入无偏位的意义上，两者都接近最佳效率。
方法一：算术（反）编码
一种实用的方法是使用 arithmetic (de)coding, as already described in the 解码您的无偏输入流。对于这个简单的案例，编写代码并不难。这是一些未优化的伪代码（咳咳，Python）执行此操作：
import random
def random_bits():
"""
Infinite generator generating a stream of random bits,
with 0 and 1 having equal probability.
"""
global bit_count # keep track of how many bits were produced
while True:
bit_count += 1
yield random.choice([0, 1])
def bernoulli(p):
"""
Infinite generator generating 1-bits with probability p
and 0-bits with probability 1 - p.
"""
bits = random_bits()
low, high = 0.0, 1.0
while True:
if high <= p:
# Generate 1, rescale to map [0, p) to [0, 1)
yield 1
low, high = low / p, high / p
elif low >= p:
# Generate 0, rescale to map [p, 1) to [0, 1)
yield 0
low, high = (low - p) / (1 - p), (high - p) / (1 - p)
else:
# Use the next random bit to halve the current interval.
mid = 0.5 * (low + high)
if next(bits):
low = mid
else:
high = mid
这是一个用法示例：
import itertools
bit_count = 0
# Generate a million deviates.
results = list(itertools.islice(bernoulli(0.1), 10**6))
print("First 50:", ''.join(map(str, results[:50])))
print("Biased bits generated:", len(results))
print("Unbiased bits used:", bit_count)
print("mean:", sum(results) / len(results))
上面给出了以下示例输出：
First 50: 00000000000001000000000110010000001000000100010000
Biased bits generated: 1000000
Unbiased bits used: 469036
mean: 0.100012
正如所承诺的那样，我们使用源无偏流中的不到五十万位生成了输出偏向流的 100 万位。
出于优化目的，将其转换为 C/C++ 时，使用基于整数的定点运算而不是浮点运算对其进行编码可能更有意义。
方法二：基于整数的算法
这里没有尝试将算术解码方法转换为直接使用整数，而是一种更简单的方法。它不再完全是算术解码，但并非完全无关，并且它实现了与上述浮点版本接近相同的输出偏置位/输入无偏置位比率。它的组织方式使所有数量都适合一个无符号的 32 位整数，因此应该很容易转换为 C/C++。该代码专门针对 p 是 1/200 的精确倍数的情况，但这种方法适用于任何可以表示为分母相当小的有理数的 p。
def bernoulli_int(p):
"""
Infinite generator generating 1-bits with probability p
and 0-bits with probability 1 - p.
p should be an integer multiple of 1/200.
"""
bits = random_bits()
# Assuming that p has a resolution of 0.05, find p / 0.05.
p_int = int(round(200*p))
value, high = 0, 1
while True:
if high < 2**31:
high = 2 * high
value = 2 * value + next(bits)
else:
# Throw out everything beyond the last multiple of 200, to
# avoid introducing a bias.
discard = high - high % 200
split = high // 200 * p_int
if value >= discard: # rarer than 1 time in 10 million
value -= discard
high -= discard
elif value >= split:
yield 0
value -= split
high = discard - split
else:
yield 1
high = split
关键的观察是，每次我们到达 while 循环的开始时，value 均匀分布在 [0, high) 中的所有整数中，并且独立于所有先前的 -输出位。如果你关心速度而不是完美的正确性，你可以去掉 discard 和 value >= discard 分支：那只是为了确保我们输出 0 和 1正确的概率。抛开这个复杂因素，您将获得 几乎 正确的概率。此外，如果您使 p 的分辨率等于 1/256 而不是 1/200，则可能耗时的除法和模运算可以替换为位运算。
使用与之前相同的测试代码，但使用 bernoulli_int 代替 bernoulli，我得到 p=0.1 的以下结果：
First 50: 00000010000000000100000000000000000000000110000100
Biased bits generated: 1000000
Unbiased bits used: 467997
mean: 0.099675
虽然这个问题已有 5 年历史，但我相信我有一些有价值的东西可以补充。虽然 SIMD 和算术解码无疑是伟大的技术，但很难忽视@mindriot 建议的 非常简单易懂。
但是，您将如何有效且快速地实施此解决方案还不是很明显。对于 256 位 (0.00390625) 的分辨率，您可以编写包含 256 个案例的 switch 语句，然后手动确定每个案例所需的布尔表达式。编写这个程序需要一段时间，但它会编译成 C/C++ 中的非常快速的跳转 table。
但是，如果您想要 2^16 位甚至 2^64 位的分辨率怎么办？后者是 5.4210109E-20 的分辨率，比我们大多数人需要的更精确。手工绝对不可能完成这项任务，但我们实际上可以构建一个小型虚拟机，仅需 30 行 C 代码即可快速完成此任务。
让我们构建 256 位分辨率的机器。我将定义 probability = resolution/256。例如，当 resolution = 64 时，则 probability = 0.25。事实证明，分子（分辨率）实际上在其二进制表示中隐式编码了所需的布尔运算。
例如，什么表达式生成 probability = 0.69140625 = 177/256？分辨率为177，二进制为10110001。让AND = 0和OR = 1。我们在第一个非零最低有效位 之后开始 并向最高有效位读取。将 0/1 映射到 AND/OR。因此，从 b1 开始，从右到左读取，我们生成布尔表达式 (((((((b1 and b2) and b3) and b4) or b5) or b6) and b7) or b8)。 computer-generated 真相 table 将确认 177 个案例产生 True。再举一个例子，probability = 0.4375 = 112/256 给出二进制的分辨率为 01110000。在第一个非零 LSB (011) 之后按顺序读取 3 位得到 ((b1 | b2) | b3) & b4).
因为我们只需要两个 AND 和 OR 操作，并且由于分辨率编码了我们需要的精确布尔表达式，因此可以对虚拟机进行编程，将分辨率解释为位码。 AND 和 OR 只是直接作用于无偏随机数生成器输出的操作码。这是我的示例 C 代码：
uint64_t rng_bias (uint64_t *state, const uint8_t resolution)
{
if (state == NULL) return 0;
//registers
uint64_t R0 = 0;
uint8_t PC = __builtin_ctz(resolution|0x80);
//opcodes
enum
{
OP_ANDI = 0,
OP_ORI = 1,
};
//execute instructions in sequence from LSB -> MSB
while (PC != (uint8_t) 0x8)
{
switch((resolution >> PC++) & (uint8_t) 0x1)
{
case OP_ANDI:
R0 &= rng_generator(state);
break;
case OP_ORI:
R0 |= rng_generator(state);
break;
}
}
return R0;
}
虚拟机无非就是2个寄存器和2个操作码。我正在使用 GCC 的内置函数 ctz 来计算尾随零位，以便我可以轻松找到第一个非零 LSB。我 bitwise-or 带有 0x80 的 ctz 参数，因为传递零是未定义的。任何其他像样的编译器都应该具有类似的功能。请注意，与我手动展示的示例不同，VM 解释从第一个非零 LSB 开始的位码，而不是 after。这是因为我需要至少调用一次 PRNG 来生成基础 p=0.5 和 p=0.0 案例。
state 指针和 rng_generator() 调用用于连接随机数生成器。例如，出于演示目的，我可以使用 Marsaglia 的 Xorshift64:
uint64_t rng_generator(uint64_t *state)
{
uint64_t x = *state;
x ^= x << 13;
x ^= x >> 7;
x ^= x << 17;
return *state = x;
}
所有 user/you 需要做的就是管理一个单独的 uint64_t state 变量，该变量必须在使用任一函数之前适当地设置种子。
缩放到 2^64 位或任何其他所需的任意分辨率非常容易。使用 ctzll 代替 unsigned long long 参数，将 uint8_t 类型更改为 uint64_t，并将 while 循环检查更改为 64 而不是 8。就是这样！现在最多调用 64 次 PRNG，速度相当快，我们可以访问 5.4210109E-20 分辨率。
这里的关键是我们实际上是免费获得位码的。没有词法分析、解析或任何其他典型的 VM 解释器任务。用户通过分辨率提供它，而没有意识到。在他们看来，只是概率的分子。对于我们这些实现者来说，它只不过是一串供我们的 VM 解释的位码。
解释位码为何有效需要一篇完全不同且更长的文章。在概率论中，问题是确定给定概率的生成事件（所有样本点的集合）。与从密度函数生成随机数的通常的逆 CDF 问题没有什么不同。从计算机科学的角度来看，在 256 位分辨率的情况下，我们正在遍历一棵深度为 8 的二叉树，其中每个节点代表一个概率。父节点是p=0.5。向左遍历表示AND次操作，向右遍历表示OR次操作。遍历和节点深度直接映射到我们之前讨论过的 LSB->MSB 位编码。
通常，随机数生成器 returns 比特流，在每个位置观察到 0 或 1 的概率是相等的（即 50%）。我们称其为无偏 PRNG。
我需要生成一串具有以下内容的伪随机位 属性：在每个位置看到 1 的概率是 p（即看到 0 的概率是 1-p）。参数p是0到1之间的实数；在我的问题中，碰巧它的分辨率为 0.5%，即它可以取值 0%、0.5%、1%、1.5%、...、99.5%、100%。
请注意，p 是概率而不是精确分数。 n位流中实际设置为1的位数必须服从二项分布B(n, p).
有一种天真的方法可以使用无偏PRNG来生成每个位的值（伪代码）：
generate_biased_stream(n, p):
result = []
for i in 1 to n:
if random_uniform(0, 1) < p:
result.append(1)
else:
result.append(0)
return result
这样的实现比生成无偏流的实现要慢得多，因为它为每个位调用一次随机数生成器函数；而无偏流生成器每个字大小调用一次（例如，它可以通过一次调用生成 32 或 64 个随机位）。
我想要一个更快的实现，即使它稍微牺牲了随机性。想到的一个想法是预先计算查找 table：对于 p 的 200 个可能值中的每一个，使用较慢的算法计算 C 8 位值并将它们保存在 table 中。然后快速算法会随机选择其中一个来生成 8 个倾斜位。
用于查看需要多少内存的粗略计算： C 至少应为 256（可能的 8 位值的数量），可能更多以避免采样效应；比方说 1024。也许这个数字应该根据 p 而变化，但让我们保持简单并假设平均值是 1024。 由于 p => 有 200 个值 => 总内存使用量为 200 KB。这还不错，可能适合 L2 缓存 (256 KB)。我还需要评估一下，看看是否有抽样效应引入偏差，这种情况下C就得调高了。
此解决方案的不足之处在于它一次只能生成 8 位，即使需要大量工作，而无偏 PRNG 只需几条算术指令即可一次生成 64 位。
我想知道是否有更快的方法，基于位操作而不是查找 tables。例如直接修改随机数生成代码，为每一位引入偏差。这将实现与无偏 PRNG 相同的性能。
编辑 3 月 5 日
谢谢大家的建议，我得到了很多有趣的想法和建议。以下是最重要的：
- 更改问题要求，使 p 的分辨率为 1/256 而不是 1/200。这允许更有效地使用位，并且还提供了更多的优化机会。我想我可以做出这个改变。
- 使用算术编码有效地消耗来自无偏生成器的位。通过上述分辨率更改，这变得容易多了。
- 一些人认为 PRNG 非常快，因此使用算术编码实际上可能会由于引入的开销而使代码变慢。相反，我应该始终使用最坏情况下的位数并优化该代码。请参阅下面的基准。
- @rici 建议使用 SIMD。这是一个不错的想法，只有当我们总是消耗固定数量的比特时它才有效。
基准（没有算术解码）
注意：正如你们许多人所建议的，我将分辨率从 1/200 更改为 1/256。
我写了 several implementations 简单地采用 8 个随机无偏位并生成 1 个有偏位的朴素方法：
- 没有 SIMD
- 使用 Agner Fog 的向量类库的 SIMD，
- 使用内部函数的 SIMD
我使用两个无偏伪随机数生成器：
- xorshift128plus
- Ranvec1 (Mersenne Twister-like) 来自 Agner Fog 的图书馆。
我还测量了无偏 PRNG 的速度以进行比较。以下是结果：
RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 16.081 16.125 16.093 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 0.778 0.783 0.812 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.176 2.184 2.145 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.129 2.151 2.183 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
与标量方法相比，SIMD 将性能提高了 3 倍。正如预期的那样，它比无偏生成器慢 8 倍。
最快的偏置生成器达到 2.1 Gb/s。
RNG: xorshift128plus
Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.300 21.486 21.483 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical : 104,857,600
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 22.660 22.661 24.662 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 1.065 1.102 1.078 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 4.972 4.971 4.970 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 4.955 4.971 4.971 [Gb/s]
Number of ones: 104,869,407 104,869,407 104,869,407
Theoretical : 104,857,600
对于 xorshift，与标量方法相比，SIMD 将性能提高了 5 倍。它比无偏生成器慢 4 倍。请注意，这是 xorshift 的标量实现。
最快的偏置生成器达到 4.9 Gb/s。
RNG: xorshift128plus_avx2
Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 18.754 21.494 21.878 [Gb/s]
Number of ones: 536,867,655 536,867,655 536,867,655
Theoretical : 104,857,600
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 54.126 54.071 54.145 [Gb/s]
Number of ones: 536,874,540 536,880,718 536,891,316
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 1.093 1.103 1.063 [Gb/s]
Number of ones: 104,868,930 104,868,930 104,868,930
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 19.567 19.578 19.555 [Gb/s]
Number of ones: 104,836,115 104,846,215 104,835,129
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 19.551 19.589 19.557 [Gb/s]
Number of ones: 104,831,396 104,837,429 104,851,100
Theoretical : 104,857,600
此实现将 AVX2 用于 运行 4 个并行的无偏异或移位生成器。
最快的偏置生成器达到 19.5 Gb/s。
算术解码基准
简单的测试表明算术解码代码是瓶颈，而不是PRNG。所以我只对最昂贵的 PRNG 进行基准测试。
RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry)
Method: Arithmetic decoding (floating point)
Gbps/s: 0.068 0.068 0.069 [Gb/s]
Number of ones: 10,235,580 10,235,580 10,235,580
Theoretical : 10,240,000
Method: Arithmetic decoding (fixed point)
Gbps/s: 0.263 0.263 0.263 [Gb/s]
Number of ones: 10,239,367 10,239,367 10,239,367
Theoretical : 10,240,000
Method: Unbiased with 1/1 efficiency (incorrect, baseline)
Gbps/s: 12.687 12.686 12.684 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical : 104,857,600
Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline)
Gbps/s: 14.536 14.536 14.536 [Gb/s]
Number of ones: 536,875,204 536,875,204 536,875,204
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency
Gbps/s: 0.754 0.754 0.754 [Gb/s]
Number of ones: 104,867,269 104,867,269 104,867,269
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=vectorclass
Gbps/s: 2.094 2.095 2.094 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
Method: Biased with 1/8 efficiency, SIMD=intrinsics
Gbps/s: 2.094 2.094 2.095 [Gb/s]
Number of ones: 104,859,067 104,859,067 104,859,067
Theoretical : 104,857,600
简单的不动点方法达到 0.25 Gb/s，而朴素的标量方法快 3 倍，朴素的 SIMD 方法快 8 倍。可能有一些方法可以进一步优化 and/or 并行化算术解码方法，但由于其复杂性，我决定就此打住并选择简单的 SIMD 实现。
谢谢大家的帮助。
如果p接近于0，则可以计算出第n位是第一个为1的位的概率；然后你计算一个介于 0 和 1 之间的随机数，并相应地选择 n。例如，如果 p = 0.005 (0.5%)，并且随机数是 0.638128，您可能会计算（我在这里猜测）n = 321，因此您填充 321 个 0 位和一个位集。
如果p接近1，用1-p代替p，设置1位加一个0位。
如果p不接近1或0，对所有256个8位序列做一个table，计算它们的累积概率，然后得到一个随机数，在数组中进行二分查找累积概率，可以设置8位。
假设 1 出现的概率是 6.25% (1/16)。 4 位数字有 16 种可能的位模式：
0000,0001, ..., 1110,1111.
现在，只需像以前一样生成一个随机数，并将半字节边界处的每个 1111 替换为 1，并将其他所有内容变为 0。
根据其他概率进行相应调整。
你可以做的一件事是多次从底层无偏生成器中采样，得到几个32位或64位的字，然后进行按位布尔运算。例如，对于 4 个单词 b1,b2,b3,b4，您可以得到以下分布：
expression | p(bit is 1)
-----------------------+-------------
b1 & b2 & b3 & b4 | 6.25%
b1 & b2 & b3 | 12.50%
b1 & b2 & (b3 | b4) | 18.75%
b1 & b2 | 25.00%
b1 & (b2 | (b3 & b4)) | 31.25%
b1 & (b2 | b3) | 37.50%
b1 & (b2 | b3 | b4)) | 43.75%
b1 | 50.00%
类似的构造可以用于更精细的分辨率。它变得有点乏味并且仍然需要更多的生成器调用，但至少不是每位一个。这类似于 a3f 的答案，但可能更容易实现，而且我怀疑比扫描 0xF nybbles 的单词更快。
请注意，对于您想要的 0.5% 的分辨率，您需要 8 个无偏见的词来对应一个有偏见的词，这将为您提供 (0.5^8) = 0.390625% 的分辨率。
给出精确结果的一种方法是首先为一个 k 位块随机生成符合二项分布的 1 位的数量，然后使用以下方法之一生成具有恰好那么多位的 k 位字方法 here。例如mic006的方法只需要大约log k k位随机数，而我的只需要一个。
它的作用
此实现通过“/dev/urandom”特殊字符文件的接口对随机设备内核模块进行单次调用，以获取表示所有值所需的随机数据数给定的决议。最大可能分辨率为 1/256^2，因此 0.005 可以表示为：
328/256^2,
即：
分辨率：256*256
x: 328
错误 0.000004883。
它是如何做到的
该实现计算位数 bits_per_byte，这是处理给定分辨率所需的均匀分布位数，即表示所有 @resolution 值。然后它会调用随机化设备（如果定义了 URANDOM_DEVICE，则调用“/dev/urandom”，否则它将通过调用“/dev/random”来使用来自设备驱动程序的额外噪声，如果存在，则可能会阻塞没有足够的比特熵）来获得所需数量的均匀分布的字节并填充随机字节数组 rnd_bytes。最后，它从 rnd_bytes 数组的每个 bytes_per_byte 字节中读取每个伯努利样本所需的位数，并将这些位的整数值与 x/resolution 给出的单个伯努利结果的成功概率进行比较。如果值命中，即它落在我们任意选择为 [0, x/resolution) 段的 x/resolution 长度段中，那么我们注意到成功并将 1 插入到结果数组中。
从随机设备读取：
/* if defined use /dev/urandom (will not block),
* if not defined use /dev/random (may block)*/
#define URANDOM_DEVICE 1
/*
* @brief Read @outlen bytes from random device
* to array @out.
*/
int
get_random_samples(char *out, size_t outlen)
{
ssize_t res;
#ifdef URANDOM_DEVICE
int fd = open("/dev/urandom", O_RDONLY);
if (fd == -1) return -1;
res = read(fd, out, outlen);
if (res < 0) {
close(fd);
return -2;
}
#else
size_t read_n;
int fd = open("/dev/random", O_RDONLY);
if (fd == -1) return -1;
read_n = 0;
while (read_n < outlen) {
res = read(fd, out + read_n, outlen - read_n);
if (res < 0) {
close(fd);
return -3;
}
read_n += res;
}
#endif /* URANDOM_DEVICE */
close(fd);
return 0;
}
填写伯努利样本向量：
/*
* @brief Draw vector of Bernoulli samples.
* @details @x and @resolution determines probability
* of success in Bernoulli distribution
* and accuracy of results: p = x/resolution.
* @param resolution: number of segments per sample of output array
* as power of 2: max resolution supported is 2^24=16777216
* @param x: determines used probability, x = [0, resolution - 1]
* @param n: number of samples in result vector
*/
int
get_bernoulli_samples(char *out, uint32_t n, uint32_t resolution, uint32_t x)
{
int res;
size_t i, j;
uint32_t bytes_per_byte, word;
unsigned char *rnd_bytes;
uint32_t uniform_byte;
uint8_t bits_per_byte;
if (out == NULL || n == 0 || resolution == 0 || x > (resolution - 1))
return -1;
bits_per_byte = log_int(resolution);
bytes_per_byte = bits_per_byte / BITS_PER_BYTE +
(bits_per_byte % BITS_PER_BYTE ? 1 : 0);
rnd_bytes = malloc(n * bytes_per_byte);
if (rnd_bytes == NULL)
return -2;
res = get_random_samples(rnd_bytes, n * bytes_per_byte);
if (res < 0)
{
free(rnd_bytes);
return -3;
}
i = 0;
while (i < n)
{
/* get Bernoulli sample */
/* read byte */
j = 0;
word = 0;
while (j < bytes_per_byte)
{
word |= (rnd_bytes[i * bytes_per_byte + j] << (BITS_PER_BYTE * j));
++j;
}
uniform_byte = word & ((1u << bits_per_byte) - 1);
/* decision */
if (uniform_byte < x)
out[i] = 1;
else
out[i] = 0;
++i;
}
free(rnd_bytes);
return 0;
}
用法：
int
main(void)
{
int res;
char c[256];
res = get_bernoulli_samples(c, sizeof(c), 256*256, 328); /* 328/(256^2) = 0.0050 */
if (res < 0) return -1;
return 0;
}
如果您准备基于 256 个可能的值来近似 p，并且您有一个 PRNG 可以生成各个位彼此独立的统一值，那么您可以使用矢量化比较从单个随机数中产生多个偏置位。
只有在 (1) 您担心随机数质量和 (2) 您可能需要大量具有相同偏差的位时才值得这样做。第二个要求似乎是由原始问题暗示的，它批评了一个提出的解决方案，如下： "A deficiency of this solution is that it can generate only 8 bits at once, even that with a lot of work, while an unbiased PRNG can generate 64 at once with just a few arithmetic instructions." 这里，暗示似乎是 useful 来生成一个一次调用中的大块偏置位。
随机数质量是一个困难的课题。即使不是不可能，也很难衡量，因此不同的人会提出不同的衡量标准，强调 and/or 贬低 "randomness" 的不同方面。通常可以权衡随机数生成的速度以获得较低的 "quality"；这是否值得做取决于您的具体应用。
最简单的随机数质量测试涉及单个值的分布和生成器的周期长度。 C 库 rand 和 Posix random 函数的标准实现通常会通过分布测试，但循环长度对于长 运行 应用程序来说是不够的。
不过，这些生成器通常非常快：random 的 glibc 实现只需要几个周期，而经典的线性同余生成器 (LCG) 需要乘法和加法。 （或者，在 glibc 实现的情况下，上面的三个生成 31 位。）如果这足以满足您的质量要求，那么尝试优化就没有意义了，尤其是在偏差概率频繁变化的情况下。
请记住，循环长度应比预期的样本数长很多；理想情况下，它应该大于该数字的平方，因此如果您希望生成千兆字节的随机数，则循环长度为 231 的线性同余生成器 (LCG) 是不合适的数据。即使是声称周期长度约为 235 的 Gnu 三项式非线性加性反馈发生器，也不应用于需要数百万样本的应用程序。
另一个更难测试的质量问题与连续样本的独立性有关。短周期长度在这个指标上完全失败，因为一旦重复开始，生成的随机数就会与历史值精确相关。 Gnu三项式算法虽然周期较长，但由于生成的ith个随机数，ith有明显的相关性=43=]ri,永远是两个值r之一i−3+ri−31 或 ri−3+ri−31+1。这可能会产生 surprising or at least puzzling 后果，尤其是对于伯努利实验。
这是一个使用 Agner Fog 的有用 vector class library, which abstracts away a lot of the annoying details in SSE intrinsics, and also helpfully comes with a fast vectorized random number generator（在 vectorclass.zip 存档中的 special.zip 中找到）的实现，它使我们能够从对 256 位 PRNG 的八次调用中生成 256 位.您可以阅读 Fog 博士关于他为何发现甚至 Mersenne twister 都存在质量问题的解释，以及他提出的解决方案；我真的没有资格发表评论，但它至少在我尝试过的伯努利实验中确实给出了预期的结果。
#include "vectorclass/vectorclass.h"
#include "vectorclass/ranvec1.h"
class BiasedBits {
public:
// Default constructor, seeded with fixed values
BiasedBits() : BiasedBits(1) {}
// Seed with a single seed; other possibilities exist.
BiasedBits(int seed) : rng(3) { rng.init(seed); }
// Generate 256 random bits, each with probability `p/256` of being 1.
Vec8ui random256(unsigned p) {
if (p >= 256) return Vec8ui{ 0xFFFFFFFF };
Vec32c output{ 0 };
Vec32c threshold{ 127 - p };
for (int i = 0; i < 8; ++i) {
output += output;
output -= Vec32c(Vec32c(rng.uniform256()) > threshold);
}
return Vec8ui(output);
}
private:
Ranvec1 rng;
};
在我的测试中，它在 260 毫秒内产生并计数了 268435456 位，或每纳秒一位。测试机是i5，所以没有AVX2； YMMV.
在实际用例中，p 有 201 个可能值，8 位阈值的计算将非常不精确。如果不需要这种不精确性，您可以调整以上内容以使用 16 位阈值，但代价是生成两倍多的随机数。
或者，您可以基于 10 位阈值手动进行矢量化，这将使您非常接近 0.5% 的增量，使用标准的位操作技巧进行矢量化阈值比较，方法是检查借用值向量与重复阈值相减的每第 10 位。与 std::mt19937_64 相结合，每个 64 位随机数平均有 6 位。
假设你可以访问一个随机位生成器，你可以生成一个值来逐位与p进行比较，并在你能证明生成的值小于时立即中止或大于或等于 p.
按照以下步骤以给定的概率在流中创建一个项目 p:
- 以二进制形式
0.开头 - 附加一个随机位；假设绘制了
1，您将得到0.1 - 如果结果（二进制表示法）可证明小于
p，则输出1 - 如果结果可证明大于或等于
p，输出一个0 - 否则（如果不能排除两者），继续步骤2。
假设p二进制表示为0.1001101...；如果此过程生成0.0、0.1000、0.10010、...中的任何一个，则该值不能再大于或等于p；如果生成0.11、0.101、0.100111、...中的任何一个，则该值不能小于p。
对我来说，看起来这种方法预期使用了大约两个随机位。对于固定 p，算术编码（如 Mark Dickinson 的回答所示）最多 每个偏置位（平均）消耗一个随机位；修改的成本 p 不清楚。
呃，伪随机数生成器一般都很快。我不确定这是什么语言（Python，也许），但是 "result.append"（几乎肯定包含内存分配）可能比 "random_uniform"（只做一点数学运算）慢.
如果你想优化这段代码的性能：
- 验证这是一个问题。优化需要一些工作，并且会使代码更难维护。除非必要，否则不要这样做。
- 剖析它。 运行 一些测试以确定代码的哪些部分实际上是最慢的。这些是您需要加速的部分。
- 进行更改，并验证它们是否确实更快。编译器非常聪明；通常清晰的代码会编译成更好的代码，比看起来更快的复杂代码。
如果您使用的是编译语言（即使是 JIT 编译语言），每次控制转移（if、while、函数调用等）都会影响性能。消除你能消除的。内存分配也（通常）非常昂贵。
如果您使用的是解释性语言，那么一切皆有可能。最简单的代码很可能是最好的。解释器的开销会让你做什么都相形见绌，所以尽可能减少它的工作。
我只能猜测你的性能问题出在哪里：
- 内存分配。以完整大小预分配数组，稍后填充条目。这确保在添加条目时不需要重新分配内存。
- 分支机构。您可以通过转换结果或类似的东西来避免 "if" 。这在很大程度上取决于编译器。检查程序集（或配置文件）以验证它是否符合您的要求。
- 数值类型。找出您的随机数生成器本机使用的类型，并在该类型中进行算术运算。例如，如果生成器自然 returns 32 位无符号整数，则先将 "p" 缩放到该范围，然后将其用于比较。
顺便说一下，如果您真的想使用尽可能少的随机数，请使用 "arithmetic coding" 来解码您的随机流。不会很快。
您将获得理论上的最佳行为，即真正地最小化使用随机数生成器并能够对任何概率建模p 完全正确， 如果你使用 arithmetic coding.
算术编码是一种数据压缩形式，它将消息表示为数字范围的子区间。它提供了理论上最优的编码，并且可以为每个输入符号使用小数位数。
这个想法是这样的：假设您有一个随机位序列，它们是 1，概率为 p。为方便起见，我将使用 q 作为该位为零的概率。 (q = 1-p)。算术编码分配给数字范围的每个位部分。对于第一位，如果输入为 0，则分配区间 [0, q)，如果输入为 1，则分配区间 [q, 1)。后续位分配当前范围的比例子区间。例如，假设 q = 1/3 输入 1 0 0 将被编码成这样：
Initially [0, 1), range = 1
After 1 [0.333, 1), range = 0.6666
After 0 [0.333, 0.5555), range = 0.2222
After 0 [0.333, 0.407407), range = 0.074074
第一个数字，1，选择范围的前三分之二（1-q）；第二个数字 0 选择 that 的后三分之一，依此类推。
在第一步和第二步之后，区间跨越中点；但是第三步之后完全在中点以下，所以可以输出第一个压缩数字：0。该过程继续，并添加一个特殊的 EOF 符号作为终止符。
这与您的问题有什么关系？压缩后的输出将具有随机的零和等概率的一。因此，要获得概率为 p 的位，只需假设你的 RNG 的输出是上述算术编码的结果，然后 对其应用解码器过程。 也就是说，读取位就好像它们将行间隔细分成越来越小的部分一样。例如，我们从 RNG 中读取 01 后，我们将在 [0.25, 0.5) 范围内。继续读取位，直到足够的输出为 "decoded"。由于您正在模拟解压缩，您将得到比输入更多的随机位。因为算术编码在理论上是最优的，所以无法将 RNG 输出转换为更多有偏差的位在不牺牲随机性的情况下：您获得了真正的最大值。
要注意的是，您无法通过几行代码完成此操作，而且我不知道有什么库可以为您指明（尽管一定有一些您可以使用）。不过，这还是很简单的。 above article provides code for a general-purpose encoder and decoder, in C. It's pretty straightforward, and it supports multiple input symbols with arbitrary probabilities; in your case a far simpler implementation is possible (as
Wikipedia 也有一个有趣的讨论算术编码被认为是基数的变化，这是查看任务的另一种方式。
从信息论的角度来看，有偏差的比特流（p != 0.5）比无偏差的比特流具有更少的信息，所以理论上它应该（平均）少于 比无偏输入的 1 位来产生有偏输出流的一位。例如，p = 0.1的伯努利随机变量的entropy是-0.1 * log2(0.1) - 0.9 * log2(0.9)位，也就是0.469位左右。这表明对于 p = 0.1 的情况，我们应该能够为每个无偏输入位生成两位多一点的输出流。
下面，我给出两种产生偏置位的方法。在需要尽可能少的输入无偏位的意义上，两者都接近最佳效率。
方法一：算术（反）编码
一种实用的方法是使用 arithmetic (de)coding, as already described in the
import random
def random_bits():
"""
Infinite generator generating a stream of random bits,
with 0 and 1 having equal probability.
"""
global bit_count # keep track of how many bits were produced
while True:
bit_count += 1
yield random.choice([0, 1])
def bernoulli(p):
"""
Infinite generator generating 1-bits with probability p
and 0-bits with probability 1 - p.
"""
bits = random_bits()
low, high = 0.0, 1.0
while True:
if high <= p:
# Generate 1, rescale to map [0, p) to [0, 1)
yield 1
low, high = low / p, high / p
elif low >= p:
# Generate 0, rescale to map [p, 1) to [0, 1)
yield 0
low, high = (low - p) / (1 - p), (high - p) / (1 - p)
else:
# Use the next random bit to halve the current interval.
mid = 0.5 * (low + high)
if next(bits):
low = mid
else:
high = mid
这是一个用法示例：
import itertools
bit_count = 0
# Generate a million deviates.
results = list(itertools.islice(bernoulli(0.1), 10**6))
print("First 50:", ''.join(map(str, results[:50])))
print("Biased bits generated:", len(results))
print("Unbiased bits used:", bit_count)
print("mean:", sum(results) / len(results))
上面给出了以下示例输出：
First 50: 00000000000001000000000110010000001000000100010000
Biased bits generated: 1000000
Unbiased bits used: 469036
mean: 0.100012
正如所承诺的那样，我们使用源无偏流中的不到五十万位生成了输出偏向流的 100 万位。
出于优化目的，将其转换为 C/C++ 时，使用基于整数的定点运算而不是浮点运算对其进行编码可能更有意义。
方法二：基于整数的算法
这里没有尝试将算术解码方法转换为直接使用整数，而是一种更简单的方法。它不再完全是算术解码，但并非完全无关，并且它实现了与上述浮点版本接近相同的输出偏置位/输入无偏置位比率。它的组织方式使所有数量都适合一个无符号的 32 位整数，因此应该很容易转换为 C/C++。该代码专门针对 p 是 1/200 的精确倍数的情况，但这种方法适用于任何可以表示为分母相当小的有理数的 p。
def bernoulli_int(p):
"""
Infinite generator generating 1-bits with probability p
and 0-bits with probability 1 - p.
p should be an integer multiple of 1/200.
"""
bits = random_bits()
# Assuming that p has a resolution of 0.05, find p / 0.05.
p_int = int(round(200*p))
value, high = 0, 1
while True:
if high < 2**31:
high = 2 * high
value = 2 * value + next(bits)
else:
# Throw out everything beyond the last multiple of 200, to
# avoid introducing a bias.
discard = high - high % 200
split = high // 200 * p_int
if value >= discard: # rarer than 1 time in 10 million
value -= discard
high -= discard
elif value >= split:
yield 0
value -= split
high = discard - split
else:
yield 1
high = split
关键的观察是，每次我们到达 while 循环的开始时，value 均匀分布在 [0, high) 中的所有整数中，并且独立于所有先前的 -输出位。如果你关心速度而不是完美的正确性，你可以去掉 discard 和 value >= discard 分支：那只是为了确保我们输出 0 和 1正确的概率。抛开这个复杂因素，您将获得 几乎 正确的概率。此外，如果您使 p 的分辨率等于 1/256 而不是 1/200，则可能耗时的除法和模运算可以替换为位运算。
使用与之前相同的测试代码，但使用 bernoulli_int 代替 bernoulli，我得到 p=0.1 的以下结果：
First 50: 00000010000000000100000000000000000000000110000100
Biased bits generated: 1000000
Unbiased bits used: 467997
mean: 0.099675
虽然这个问题已有 5 年历史，但我相信我有一些有价值的东西可以补充。虽然 SIMD 和算术解码无疑是伟大的技术，但很难忽视@mindriot 建议的
但是，您将如何有效且快速地实施此解决方案还不是很明显。对于 256 位 (0.00390625) 的分辨率，您可以编写包含 256 个案例的 switch 语句，然后手动确定每个案例所需的布尔表达式。编写这个程序需要一段时间，但它会编译成 C/C++ 中的非常快速的跳转 table。
但是，如果您想要 2^16 位甚至 2^64 位的分辨率怎么办？后者是 5.4210109E-20 的分辨率，比我们大多数人需要的更精确。手工绝对不可能完成这项任务，但我们实际上可以构建一个小型虚拟机，仅需 30 行 C 代码即可快速完成此任务。
让我们构建 256 位分辨率的机器。我将定义 probability = resolution/256。例如，当 resolution = 64 时，则 probability = 0.25。事实证明，分子（分辨率）实际上在其二进制表示中隐式编码了所需的布尔运算。
例如，什么表达式生成 probability = 0.69140625 = 177/256？分辨率为177，二进制为10110001。让AND = 0和OR = 1。我们在第一个非零最低有效位 之后开始 并向最高有效位读取。将 0/1 映射到 AND/OR。因此，从 b1 开始，从右到左读取，我们生成布尔表达式 (((((((b1 and b2) and b3) and b4) or b5) or b6) and b7) or b8)。 computer-generated 真相 table 将确认 177 个案例产生 True。再举一个例子，probability = 0.4375 = 112/256 给出二进制的分辨率为 01110000。在第一个非零 LSB (011) 之后按顺序读取 3 位得到 ((b1 | b2) | b3) & b4).
因为我们只需要两个 AND 和 OR 操作，并且由于分辨率编码了我们需要的精确布尔表达式，因此可以对虚拟机进行编程，将分辨率解释为位码。 AND 和 OR 只是直接作用于无偏随机数生成器输出的操作码。这是我的示例 C 代码：
uint64_t rng_bias (uint64_t *state, const uint8_t resolution)
{
if (state == NULL) return 0;
//registers
uint64_t R0 = 0;
uint8_t PC = __builtin_ctz(resolution|0x80);
//opcodes
enum
{
OP_ANDI = 0,
OP_ORI = 1,
};
//execute instructions in sequence from LSB -> MSB
while (PC != (uint8_t) 0x8)
{
switch((resolution >> PC++) & (uint8_t) 0x1)
{
case OP_ANDI:
R0 &= rng_generator(state);
break;
case OP_ORI:
R0 |= rng_generator(state);
break;
}
}
return R0;
}
虚拟机无非就是2个寄存器和2个操作码。我正在使用 GCC 的内置函数 ctz 来计算尾随零位，以便我可以轻松找到第一个非零 LSB。我 bitwise-or 带有 0x80 的 ctz 参数，因为传递零是未定义的。任何其他像样的编译器都应该具有类似的功能。请注意，与我手动展示的示例不同，VM 解释从第一个非零 LSB 开始的位码，而不是 after。这是因为我需要至少调用一次 PRNG 来生成基础 p=0.5 和 p=0.0 案例。
state 指针和 rng_generator() 调用用于连接随机数生成器。例如，出于演示目的，我可以使用 Marsaglia 的 Xorshift64:
uint64_t rng_generator(uint64_t *state)
{
uint64_t x = *state;
x ^= x << 13;
x ^= x >> 7;
x ^= x << 17;
return *state = x;
}
所有 user/you 需要做的就是管理一个单独的 uint64_t state 变量，该变量必须在使用任一函数之前适当地设置种子。
缩放到 2^64 位或任何其他所需的任意分辨率非常容易。使用 ctzll 代替 unsigned long long 参数，将 uint8_t 类型更改为 uint64_t，并将 while 循环检查更改为 64 而不是 8。就是这样！现在最多调用 64 次 PRNG，速度相当快，我们可以访问 5.4210109E-20 分辨率。
这里的关键是我们实际上是免费获得位码的。没有词法分析、解析或任何其他典型的 VM 解释器任务。用户通过分辨率提供它，而没有意识到。在他们看来，只是概率的分子。对于我们这些实现者来说，它只不过是一串供我们的 VM 解释的位码。
解释位码为何有效需要一篇完全不同且更长的文章。在概率论中，问题是确定给定概率的生成事件（所有样本点的集合）。与从密度函数生成随机数的通常的逆 CDF 问题没有什么不同。从计算机科学的角度来看，在 256 位分辨率的情况下，我们正在遍历一棵深度为 8 的二叉树，其中每个节点代表一个概率。父节点是p=0.5。向左遍历表示AND次操作，向右遍历表示OR次操作。遍历和节点深度直接映射到我们之前讨论过的 LSB->MSB 位编码。