我非常接近@jwimberley 的想法，并认为我会在这里分享我的结果。我创建了一个执行以下操作的类：
- 构造函数参数：
-
CDF（标准化或非标准化），即
整合 PDF。
- 分布的下限和上限
- （可选）指示我们应该采用多少 CDF 采样点的分辨率。
- 从 CDF 计算映射 -> 随机数 x。这是我们的逆 CDF 函数。
- 通过以下方式生成随机点：
- 使用
std::random 在(0, 1] 之间生成随机概率p。
- 在我们的映射中二进制搜索对应于 p 的 CDF 值。返回与 CDF 一起计算的 x。提供了附近“桶”之间的可选线性积分，否则我们将得到 n == 分辨率离散步数。
代码：
// sampled_distribution.hh
#ifndef SAMPLED_DISTRIBUTION
#define SAMPLED_DISTRIBUTION
#include <algorithm>
#include <vector>
#include <random>
#include <stdexcept>
template <typename T = double, bool Interpolate = true>
class Sampled_distribution
{
public:
using CDFFunc = T (*)(T);
Sampled_distribution(CDFFunc cdfFunc, T low, T high, unsigned resolution = 200)
: mLow(low), mHigh(high), mRes(resolution), mDist(0.0, 1.0)
{
if (mLow >= mHigh) throw InvalidBounds();
mSampledCDF.resize(mRes + 1);
const T cdfLow = cdfFunc(low);
const T cdfHigh = cdfFunc(high);
T last_p = 0;
for (unsigned i = 0; i < mSampledCDF.size(); ++i) {
const T x = i/mRes*(mHigh - mLow) + mLow;
const T p = (cdfFunc(x) - cdfLow)/(cdfHigh - cdfLow); // normalising
if (! (p >= last_p)) throw CDFNotMonotonic();
mSampledCDF[i] = Sample{p, x};
last_p = p;
}
}
template <typename Generator>
T operator()(Generator& g)
{
T cdf = mDist(g);
auto s = std::upper_bound(mSampledCDF.begin(), mSampledCDF.end(), cdf);
auto bs = s - 1;
if (Interpolate && bs >= mSampledCDF.begin()) {
const T r = (cdf - bs->prob)/(s->prob - bs->prob);
return r*bs->value + (1 - r)*s->value;
}
return s->value;
}
private:
struct InvalidBounds : public std::runtime_error { InvalidBounds() : std::runtime_error("") {} };
struct CDFNotMonotonic : public std::runtime_error { CDFNotMonotonic() : std::runtime_error("") {} };
const T mLow, mHigh;
const double mRes;
struct Sample {
T prob, value;
friend bool operator<(T p, const Sample& s) { return p < s.prob; }
};
std::vector<Sample> mSampledCDF;
std::uniform_real_distribution<> mDist;
};
#endif
以下是一些结果图。对于给定的PDF，我们需要先通过积分来解析计算CDF。
对数线性
正弦曲线
您可以通过以下演示自己尝试一下：
// main.cc
#include "sampled_distribution.hh"
#include <iostream>
#include <fstream>
int main()
{
auto logFunc = [](double x) {
const double k = -1.0;
const double m = 10;
return x*(k*std::log(x) + m - k); // PDF(x) = k*log(x) + m
};
auto sinFunc = [](double x) { return x + std::cos(x); }; // PDF(x) = 1 - sin(x)
std::mt19937 gen;
//Sampled_distribution<> dist(logFunc, 1.0, 1e4);
Sampled_distribution<> dist(sinFunc, 0.0, 6.28);
std::ofstream file("d.txt");
for (int i = 0; i < 100000; i++) file << dist(gen) << std::endl;
}
数据是用python绘制的。
// dist_plot.py
import numpy as np
import matplotlib.pyplot as plt
d = np.loadtxt("d.txt")
fig, ax = plt.subplots()
bins = np.arange(d.min(), d.max(), (d.max() - d.min())/50)
ax.hist(d, edgecolor='white', bins=bins)
plt.show()
运行演示：
clang++ -std=c++11 -stdlib=libc++ main.cc -o main; ./main; python dist_plot.py