Wie erstellt man eine benutzerdefinierte Zufallsverteilung?

Question

Normalerweise erzeugen I-Werte mit the built in random functions, aber jetzt brauche ich eine zufällige Verteilung der Form zu schaffen

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

Ist es möglich, eine benutzerdefinierte Zufallsverteilungsfunktion zu definieren? Für mein aktuelles Modell habe ich x = [1, 1.4e7), k = -0.905787102751, m = 14.913170454. Im Idealfall würde Ich mag es zu arbeiten, wie die aktuellen Einbau-Verteilungen zu tun:

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

aussehen Ich folgte Idee @ jwimberley so ziemlich auf den Punkt, und dachte, ich würde hier meine Ergebnisse teilen. Ich habe eine Klasse, führt Folgendes aus:

1. Konstruktorargumente:
   • CDF (normiert oder nicht-normalisierten), das ist die integral des PDF.
   • Untere und obere Grenze der Verteilung
   • (optional) Auflösung, die angibt, wie viele Beispielpunkte der CDF wir nehmen sollten.
2. Berechnen Sie eine Zuordnung von CDF -> Zufallszahl x. Dies ist unsere inverse CDF-Funktion.
3. Generieren Sie einen beliebigen Punkt durch:
   • eine zufällige Wahrscheinlichkeit generieren p zwischen (0, 1]std::random verwenden.
   • Binärsuche in unserem Mapping für CDF-Wert, der p entspricht. Gebe die x zurück, die zusammen mit der CDF berechnet wurde. Optionale lineare Integration zwischen benachbarten "Buckets" wird bereitgestellt, andernfalls erhalten wir n == Auflösung diskrete Schritte.

Der Code:

// 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 

Hier sind einige Grundstücke der Ergebnisse. Für eine gegebene PDF müssen wir zuerst die CDF durch Integration analytisch berechnen.

Log-lineare

Sinusförmige

Sie können dies ausprobieren, sich mit der folgenden Demo:

// 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; 
} 

Die Daten werden mit Python aufgetragen.

// 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() 

Führen Sie die Demo mit:

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

Answer 2

Das ist sehr gut möglich, aber es ist so viel von einem mathematischen Problem als C++ Problem. Die allgemeinste Methode zum Erzeugen eines Pseudozufallszahlengenerators ist Inverse transform sampling. Im Wesentlichen ist die CDF eines PDFs einheitlich zwischen 0 und 1 verteilt (wenn dies nicht offensichtlich ist, denken Sie daran, dass der Wert der CDF eine Wahrscheinlichkeit ist und denken Sie darüber nach). Sie müssen also einfach eine zufällige Einheitszahl zwischen 0 und 1 abtasten und die Umkehrung der CDF anwenden.

In Ihrem Fall, mit $ f (x) = k * log (x) + m $ (Sie haben keine Grenzen angegeben, aber ich nehme an, sie sind zwischen 1 und einige positive Zahl> 1) die CDF und seine Inverse sind ziemlich unordentlich - ein Problem, das ich Ihnen überlassen! Die Implementierung in C++ wird wie folgt aussehen

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

die Erzeugung Code Dann etwas wie

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

Wie an anderer Stelle erwähnt, ist das Standardverfahren für jede PDF Abtasten seiner CDF an einem Punkt aus dem Intervall gleichmäßig zufällig ausgewählt invertieren [0,1] .

Im Falle Ihres speziellen Problems ist die CDF eine einfache Funktion, aber ihre Umkehrung ist nicht. In diesem Fall kann man es unter Verwendung herkömmlicher numerischer Werkzeuge, wie Newton-Raphson-Iteration, invertieren. Leider haben Sie den Bereich für x oder die zulässigen Optionen für die Parameter m und k nicht angegeben. Ich habe dies für allgemeine m, k implementiert, und Bereich (and posted it on code review), um die C++ RandomNumberDistribution concept zu erfüllen.