logo标签列表

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

c++c++11mathrandom
4

通常我使用内置的随机函数生成值，但现在我需要创建一个形式为

的随机分布

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);
}
1
这个问题涉及到的数学和C++代码一样重要。例如，请参阅https://en.wikipedia.org/wiki/Inverse_transform_sampling。 - jwimberley
1
域名是什么？ - user1196549
@YvesDaoust 对于最初的问题，它是在1到1.4e7之间。我已经添加了一个解决方案的答案。 - pingul
请指定参数m和k的预期范围以及范围。特别地，是否考虑x小于1？ - Walter
@Walter 我已经在问题中添加了我的实际模型值。谢谢。 - pingul
4个回答
6
这是完全可能的，但它既是一个数学问题又是一个C++问题。创建伪随机数生成器最常见的方法是反函数变换采样。基本上，任何PDF的CDF均匀分布在0到1之间（如果不明确，请记住CDF的值是概率并仔细思考）。因此，你只需要在0到1之间采样一个随机均匀数，并应用反函数变换。

对于$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);
     }
}
这是非常好的建议，让我走上了正确的道路。已点赞。我添加了一个答案，概述了我如何实现它 - 这是你想要的吗？如果您觉得有任何问题，请提出改进意见。 - pingul
6

我基本上完全按照 @jwimberley 的想法，并想在这里分享我的结果。 我创建了一个类，它执行以下操作：

  1. 构造函数参数：
    • CDF（归一化或非归一化），它是PDF的积分。
    • 分布的下限和上限
    • （可选）决定我们应该取多少CDF样本点的分辨率。
  2. 计算从CDF -> 随机数 x的映射。这是我们的反向CDF函数。
  3. 通过以下方式生成随机点：
    • 使用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

以下是翻译的结果：

这里是一些结果的图表。对于给定的概率密度函数，我们需要首先通过积分来解析计算累积分布函数。

对数线性 对数线性分布

正弦波形 正弦波形分布

您可以使用以下演示来尝试：

// 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
1
关于这段代码，有很多话可以说，但这确实属于代码审查。 - Walter
1
@Walter 这篇帖子并没有要求评论。这是我回答自己问题的方式，解释了如何创建一个自定义随机分布。我对被踩感到非常惊讶。 - pingul
1
@Walter，无需测试单调性，有效的 CDF始终是单调不降的，因为任何有效的PDF必须是非负的。但对于离散分布，您必须修改二分搜索。 - pjs
我之前实现过类似的东西。但我做了一些不同的事情。首先，我的类接受一个 lambda 函数来计算 PDF 并数值地计算 CDF 或可选的 lambda 函数 CDF。这对于许多代码的客户来说更加实际，因为他们可能没有解析的 CDF。其次，我同意二分查找不是最有效的反演方法。我使用了类似于牛顿-拉弗森算法的方法。数值积分和牛顿-拉弗森都在 GSL 库中实现，以避免手动操作。 - jwimberley
@Walter，我很感激pingul所做的努力，展示了一个完整的解决方案来回答这个问题。如果你不喜欢代码风格，你可以自由地建议具体改进什么，甚至更好的方法是自己编辑答案。 - HAL9000
显示剩余4条评论
2

我非常喜欢这里介绍的概念，它们导致了一个非常精简但功能强大的生成器。我刚刚做了一些清理工作，嵌入了C ++17功能，本来想编辑pingul的答案，但结果有很大不同，所以我将其单独发布。

#pragma once

#include <algorithm>
#include <vector>
#include <random>
#include <stdexcept>

template <typename T = double, bool Interpolate = true>
class SampledDistribution {
  struct Sample { 
    T prob, value; 
    Sample(const T p, const T v): prob(p), value(v) {}
    friend bool operator<(T p, const Sample& s) { return p < s.prob; }
  };

  std::vector<Sample> SampledCDF;

public:
  struct InvalidBounds:   std::runtime_error { using std::runtime_error::runtime_error; };
  struct CDFNotMonotonic: std::runtime_error { using std::runtime_error::runtime_error; };

  template <typename F>
  SampledDistribution(F&& cdfFunc, const T low, const T high, const unsigned resolution = 256) {
    if (low >= high) throw InvalidBounds("");
    SampledCDF.reserve( resolution );
    const T cdfLow = cdfFunc(low);
    const T cdfHigh = cdfFunc(high);
    for (unsigned i = 0; i < resolution; ++i) {
      const T x = (high - low)*i/(resolution-1) + low;
      const T p = (cdfFunc(x) - cdfLow)/(cdfHigh - cdfLow); // normalising 
      if (p < SampledCDF.back()) throw CDFNotMonotonic("");
      SampledCDF.emplace_back(p, x);
    }
  }

  template <typename Engine>
  T operator()(Engine& g) {
    const T cdf = std::uniform_real_distribution<T>{0.,1.}(g);
    auto s = std::upper_bound(SampledCDF.begin(), SampledCDF.end(), cdf);
    if (Interpolate && s != SampledCDF.begin()) { 
      auto bs = s - 1;
      const T r = (cdf - bs->prob)/(s->prob - bs->prob);
      return r*bs->value + (1 - r)*s->value;
    }
    return s->value;
  }
};

这里是一个测试主函数：

#include <iostream>
#include "SampledDistribution.hpp"

int main() {
  std::mt19937 gen;
  auto sinFunc = [](double x) { return x + std::cos(x); }; // PDF(x) = 1 - sin(x)

  unsigned resolution = 32;
  std::vector<int> v(resolution,0);
  SampledDistribution dist(sinFunc, 0.0, 6.28, resolution);

  for (int i = 0; i < 100000; i++) 
    ++v[ static_cast<size_t>(dist(gen)/(6.28) * resolution) ];

  for (auto i: v)
    std::cout << i << '\t' << std::string(i/100, '*') << std::endl;

  return 0;
}

样例输出：

$ g++ -std=c++17 main.cpp && ./a.out
2882    ****************************
2217    **********************
1725    *****************
1134    ***********
690     ******
410     ****
182     *
37  
34  
162     *
411     ****
753     *******
1163    ***********
1649    ****************
2157    *********************
2796    ***************************
3426    **********************************
4048    ****************************************
4643    **********************************************
5193    ***************************************************
5390    *****************************************************
5796    *********************************************************
5979    ***********************************************************
6268    **************************************************************
6251    **************************************************************
6086    ************************************************************
5783    *********************************************************
5580    *******************************************************
5111    ***************************************************
4646    **********************************************
3964    ***************************************
3434    **********************************
0

正如其他地方指出的那样，对于任何PDF的抽样标准方法是在区间[0,1]中均匀随机选择一个点并反转其CDF。

对于您特定的问题，CDF是一个简单的函数，但它的反函数不是。在这种情况下，可以使用传统的数值工具（例如牛顿-拉弗森迭代）来反转它。不幸的是，您未指定x的范围或参数m和k的允许选择。我已经实现了通用的m、k和范围（并将其发布在代码审查中），以满足C++ RandomNumberDistribution概念。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接