概率规划

Question

在我们需要生成概率的情况下，例如具有 75% 的抛头和 25% 的抛尾的偏差硬币。按照惯例，我会这样做：

#include <cstdlib>
#include <iostream>
#include <ctime>
using namespace std;

int main()
{
    int heads=0, tails=0;
    srand(time(NULL));
    number = rand() % 100 + 1;  //Generate random number 1 to 100
          if (number <= 75) //75% chance
                heads++; //This is head
          else
                tails++; //This is tail
}

这是一个工作代码，但是当我在 SO 中为另一个用户回答类似的关于偏差硬币的问题时，一些用户提到了 100 的倍数。由于随机函数产生均匀分布，我觉得上面的代码足以模拟概率事件。

在过去的 SO 帖子中，用户 Bathsheba 提到了一些关于 100 的倍数的问题：Program that simulates a coin toss with a bias coin 我想知道与此相关的代码中可能存在哪些问题。

我的问题是：上面的代码是可以接受的代码来创建概率模拟吗？或者这些代码中是否有任何缺陷会影响模拟结果的准确性。如果上述代码存在缺陷，那么实现概率模拟的正确方法是什么？

编辑：进行了 10,000,000 次抛掷的模拟测试。它总是以大约 75.01%-75.07% 的概率产生抛头。那么当它产生一个看似准确的结果时会出现什么问题。 （生成的结果似乎没有偏差）

Answer 1

上面的代码是一个可接受的代码来创建一个模拟 可能性？或者这些代码中是否有任何缺陷会影响 模拟结果的准确性？

如果这是“可接受的”，这取决于您对可接受的定义。这肯定是不正确的，因为 operator % 会使您的概率出现偏差，因为 RAND_MAX 是 rand() 的最大值不能等于 k ​​ 100 + 99 这导致如果您想象 0-RAND_MAX 字符串的 100 长度部分那么你可以看到最后一部分可能不会产生一个完整的范围 0-99，所以你有更多的数字可以产生 0, 1, 2..., x 但不是必需的 x + 1, ..., 98, 99（0, 1, 2, ..., x 中的每个数字再出现 1 次）。这种方法的不准确性随着不均匀划分范围的更大除数而增加。

如果上面的代码有缺陷，正确的方法是什么 实现概率模拟？

您可以使用 boost，或者如果您可以运行 C++11，那么您可以使用标准库的 uniform_int_distribution。

Answer 2

由于数字的有限性质，您总是会得到有偏差的结果（增加随机数生成器的结果数量会提高准确性）

在您的样本中，您可能对什么是 75% 有更好的定义：

int main()
{
    int heads=0, tails=0;
    srand(time(NULL));
    const std::size_t Samples = 10000000;
    for(std::size_t i = 0; i < Samples; ++i) {
        int head_limit = RAND_MAX * 0.75;
        int number = rand();
        if (number <= head_limit) heads++;
        else tails++;
    }
    // heads: 7498728 [0.749873%]
    // tails: 2501272 [0.250127%]
    std::cout 
        << "heads: " << heads << " [" << double(heads) / Samples << "%]\n"
        << "tails: " << tails << " [" << double(tails) / Samples << "%]\n";
}

Answer 3

使用rand() % 100 + 1 不能像“同时生成 100 个随机数 - 恰好 75 个数字将小于 75”这样的方式工作

以其他方式 - 它不保证在 100 个随机生成的数字中，75 个数字将小于 75！

Answer 4

上面的代码是一个可接受的代码来创建一个模拟 可能性？或者这些代码中是否有任何缺陷会影响 模拟结果的准确性。

我不知道你对“可接受”的定义。但是，我会完全避免使用 rand()，例如参见 rand() Considered Harmful。

如果上述代码有缺陷，那么实现模拟的正确方法是什么？ 概率？

我会使用 std::bernoulli_distribution 和 Mersenne Twister engine。它质量高、速度快（根据 Stephan T. Lavavej 的介绍）和标准。

顺便说一下，std::bernoulli_distribution 的示例代码给出了 1/4 的时间“真”和 3/4 的“假”。 ;)

Answer 5

std::rand()generates a number between 0 and RAND_MAX，即guaranteed to be at least 32767。

假设RAND_MAX 定义为32767，rand()%100 会产生平坦分布吗？不会。从 0 到 32699，从 0 到 99 的每个值将出现 327 次。但是从 32700 到 32767，值 0 到 67 出现了一次，而 68-99 出现了 0 次。所以你的分布有 328 个 00-67 和 327 个 68-99。

此外，除非您以某种方式指定 RAND_MAX 是什么，或者在代码中使用它，否则您将受制于编译器实现对 RAND_MAX 使用的任何东西，并且您的分布将以某种未知的方式出现偏差。

如果您希望硬币在四分之三的时间内出现正面，请考虑以下情况：

if((double)std::rand()/3.0 > (double)RAND_MAX/4.0)

（如果 a > 3/4 * b 则 a/3 > b/4）。这几乎是公平的； RAND_MAX 不太可能整齐地分成 4。但它会比原始代码中的 1/327 偏差更好。

但更好的是，使用更好的随机数生成器来设置限制。

Answer 6

为了确保 1 到 100 或 0 到 99 之间的每个数字具有 P=1/100 的概率，以确保您有准确的概率排序，

然后不是使用随机生成的数字，而是使用1000个1-100的列表均匀分布，然后每次需要重新使用它们时，使用相同的随机数生成器将它们洗牌，

所以首先我们建立列表：

const int SIZE = 1000;
srand(time(NULL));
int randList [SIZE];

然后我们填充它：

void init (int randList[], const int SIZE)
{
    for (int i=0; i<SIZE; i++)
        randList[i] = i % 100;
}

然后在每 1000 次硬币试验之前，我们将列表洗牌：

void shuffle (int randList[], const int SIZE)
{
    for (int i=0; i<SIZE; i++)
        swap(randList,i,(rand() % SIZE));
}

void swap (int randList[], int a, int b)
{
    int t = randList[a];
    randList[a] = randList[b];
    randList[b] = t;
}

那么我们可以像这样进行试验：

bool trial (int randList[], const int SIZE, int trialCount)
{
    return (randList[trialCount % SIZE] < 75); // Head = True = 75%
}

然后是一组试验：

void test (bool * resultList, const int resultSize)
{
    const int SIZE = 1000;
    srand(time(NULL));
    int randList [SIZE];

    init(randList,SIZE);

    for (int i=0; i<resultSize; i++)
    {
        if (i%SIZE == 0)
            shuffle(randList,SIZE);

        resultList[i] = trial(randList,SIZE,i);
    }
}

最后，在 main 函数中我们直接使用 test 函数：

int main ()
{
    const int resultSize = 2000000; // 2 Million

    bool * resultList = new bool[resultSize];

    test(resultList,resultSize);

    // check sequence of outcomes

    return 0;
}