如何创建自定义随机分布函数？

Question

通常我使用内置的随机函数生成值，但现在我需要创建以下形式的随机分布

f(x) = k*log(x) + m

是否可以定义自定义随机分布函数？对于我的实际模型，我有x = [1, 1.4e7), k = -0.905787102751, m = 14.913170454. 理想情况下，我希望它能够像当前内置发行版一样工作：

int main() 
{
    std::mt19937 generator;

    std::uniform_real_distribution<> dist(0.0, 1.0);
    my_distribution my_dist(0.0, 10.0); // Distribution using f(x)

    double uni_val = dist(generator);
    double log_val = my_dist(generator);
}

Answer 1

这是很有可能的，但它既是一个数学问题，又是一个 C++ 问题。创建伪随机数生成器的最通用方法是逆变换采样。本质上，任何 PDF 的 CDF 均匀分布在 0 和 1 之间（如果这不明显，只需记住 CDF 的值是一个概率并思考这一点）。因此，您只需采样 0 到 1 之间的随机均匀数并应用 CDF 的逆。

在你的情况下，使用 $f(x) = k*log(x) + m$ （你没有指定界限，但我假设它们在 1 和某个正数 > 1 之间）CDF 及其逆非常混乱- 我把这个问题留给你了！C++ 中的实现如下所示

double inverseCDF(double p, double k, double m, double lowerBound, double upperBound) {
     // do math, which might include numerically finds roots of equations
}

那么生成代码将类似于

class my_distribution {
     // ... constructor, private variables, etc.
     template< class Generator >
     double operator()( Generator& g ) {
          std::uniform_real_distribution<> dist(0.0, 1.0);
          double cdf = dist(g);
          return inverseCDF(cdf,this->k,this->m,this->lowerBound,this->upperBound);
     }
}

Answer 2

我完全遵循@jwimberley 的想法，并想在这里分享我的结果。我创建了一个执行以下操作的类：

1. 构造函数参数：
   • CDF（标准化或非标准化），它是 PDF的积分。
   • 分布的下限和上限
   • （可选）指示我们应该采取多少个 CDF 样本点的分辨率。
2. 计算 CDF -> 随机数x的映射。这是我们的逆 CDF 函数。
3. 通过以下方式生成随机点：
   • 使用之间生成随机概率p。(0, 1]std::random
   • 在我们的映射中二分搜索对应于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