C++标准库中std::poisson_distribution存在Bug？

Question

我想我遇到了C++标准库中std::poisson_distribution函数的错误行为。

问题:

1. 您能确认这确实是一个bug而不是我的错误吗？
2. 假设它确实是一个bug，那么标准库中poisson_distribution函数的代码有什么问题？

细节:

以下C++代码（文件poisson_test.cc）用于生成泊松分布的数字:

#include <array>
#include <cmath>
#include <iostream>
#include <random>

int main() {
  // The problem turned out to be independent on the engine
  std::mt19937_64 engine;

  // Set fixed seed for easy reproducibility
  // The problem turned out to be independent on seed
  engine.seed(1);
  std::poisson_distribution<int> distribution(157.17);

  for (int i = 0; i < 1E8; i++) {
    const int number = distribution(engine);
    std::cout << number << std::endl;
  }
}

我按以下方式编译这段代码：

clang++ -o poisson_test -std=c++11 poisson_test.cc
./poisson_test > mypoisson.txt

以下是用于分析文件“mypoisson.txt”中随机数序列的Python脚本：

import numpy as np
import matplotlib.pyplot as plt

def expectation(x, m):
    " Poisson pdf " 
    # Use Ramanujan formula to get ln n!
    lnx = x * np.log(x) - x + 1./6. * np.log(x * (1 + 4*x*(1+2*x))) + 1./2. * np.log(np.pi)
    return np.exp(x*np.log(m) - m - lnx)

data = np.loadtxt('mypoisson.txt', dtype = 'int')

unique, counts = np.unique(data, return_counts = True)
hist = counts.astype(float) / counts.sum()
stat_err = np.sqrt(counts) / counts.sum()
plt.errorbar(unique, hist, yerr = stat_err, fmt = '.', \
             label = 'Poisson generated \n by std::poisson_distribution')
plt.plot(unique, expectation(unique, expected_mean), \
         label = 'expected probability \n density function')
plt.legend()
plt.show()

# Determine bins with statistical significance of deviation larger than 3 sigma
deviation_in_sigma = (hist - expectation(unique, expected_mean)) / stat_err
d = dict((k, v) for k, v in zip(unique, deviation_in_sigma) if np.abs(v) > 3.0)
print d

脚本生成以下图表：

你可以用肉眼看出问题。n = 158处的偏差在统计学上是显著的，实际上是22σ的偏差！

前一个图表的局部放大。

- Dmytro Oliinychenko

---

是哪个标准库呢？据我所知，Clang 在 Linux 上倾向于使用 libstdc++ ，在 Mac 上则倾向于使用 libc++。 - chris

@chris 我正在使用Ubuntu，我已经检查过它是libstdc++。打印出____GLIBCXX____的结果是20160609。关于clang版本，“clang -v”显示“clang version 3.8.0-2ubuntu4 (tags/RELEASE_380/final)”。你能在Mac上重现这个错误吗？ - Dmytro Oliinychenko

我已经使用Visual C++ 2017（32位版本）编译并运行，但我没有得到异常值。（158的值约为0.03166，略低于157的值0.03182） - 1201ProgramAlarm

仅对数学进行评论，泊松分布不是连续的而是离散的，因此没有概率密度函数。但是你可以计算概率质量函数。我知道你的线图有助于引导眼睛，但是泊松分布确实具有非整数变量，因此绘制连续线条有点误导性。 - Jan Christoph Terasa

看源代码：// NB: This case not in the book, nor in the Errata, but should be ok... - 除了一些非常基本的大学关于接受/拒绝算法的东西，我对手头的问题一无所知，但这种声明让我感到紧张... :o) - Matteo Italia

Answer 1

我的系统如下所示（Debian testing）:

libstdc++-7-dev:
  Installed: 7.2.0-16

libc++-dev:
  Installed: 3.5-2

clang:
  Installed: 1:3.8-37

g++:
  Installed: 4:7.2.0-1d1

当使用libstdc++时，我可以确认这个bug:

g++ -o pois_gcc -std=c++11 pois.cpp
clang++ -o pois_clang -std=c++11 -stdlib=libstdc++ pois.cpp
clang++ -o pois_clang_libc -std=c++11 -stdlib=libc++ pois.cpp

结果：

- Jan Christoph Terasa

---

你是否报告了那个漏洞？如果还没有的话，你应该报告一下！然后在你的问题中提供它的网址！ - Basile Starynkevitch

这个功劳不属于我。如果 @DmytroOliinychenko 想报告，我会等他。尽管如果他决定不这样做，我也可以这么做。 - Jan Christoph Terasa

已向GCC错误跟踪器报告。 - Jan Christoph Terasa

@ChristophTerasa请链接到错误报告。 - M.M

GCC漏洞跟踪链接：https://gcc.gnu.org/bugzilla/show_bug.cgi?id=83237 - Jan Christoph Terasa

@ChristophTerasa 谢谢！这个答案非常有帮助，因为现在我知道只需使用“-stdlib=libc++”编译即可获得正确的分发。我很想将其作为最终正确答案接受，但作为一名科学家，我想知道libstdc++实现有什么问题。源代码中的那条注释（// NB：此情况不在书中，也不在勘误表中，但应该没问题...）确实令人怀疑，但我无法证明它是导致错误的原因。 - Dmytro Oliinychenko