カスタムランダム分布関数の作成方法は？

Question

は、通常、私はthe built in random functionsを使用して値を生成するが、今は

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

のようになります。私はポイントにほとんどの@ jwimberleyの考え方を踏襲し、私はここに私の結果を共有するだろうと思いました。

1. コンストラクタ引数：
   • CDF（正規化または非正規化）、PDFの 積分である私は、次のないクラスを作成しました。
   • ディストリビューションの上限と下限
   • （オプション）CDFのサンプルポイント数を示す解像度。
2. CDF→乱数xからのマッピングを計算します。これは逆CDF関数です。
3. によってランダム点を生成します。
   • ランダム確率Pstd::randomを使用して(0, 1]間を生成します。
   • 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() 

実行とデモ：ASが他の場所で指摘

clang++ -std=c++11 -stdlib=libc++ main.cc -o main; ./main; python dist_plot.py 

Answer 2

これは非常に可能ですが、C++の問題として、数学の問題のその限り。疑似乱数ジェネレータを作成する最も一般的な方法はInverse transform samplingです。基本的に、PDFのCDFは0と1の間で均等に分布しています（明らかでない場合は、CDFの値が確率であり、これを考えることを忘れないでください）。したがって、0と1の間のランダムな一様数をサンプリングし、CDFの逆数を適用するだけで済みます。

$ f（x）= k * log（x）+ m $（境界は指定していませんが、1と正の数> 1の間であると仮定します）逆はかなり面倒です - 私はあなたに残す問題です！ 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 3

、任意のPDFをサンプリングするための標準的な方法は、間隔から一様にランダム選択された時点で、そのCDFを反転することである[0,1] 。

特に問題がある場合、CDFは単純な関数ですが、逆関数はありません。この場合、Newton-Raphson反復などの従来の数値ツールを使用して逆転させることができます。残念ながら、xの範囲、またはパラメータmとkの許容される選択肢の指定に失敗しました。私はこれを一般的なm,k、および範囲（and posted it on code review）に実装して、C++ RandomNumberDistribution conceptを満足させました。