使用所谓的接受拒绝方法可以最轻松地完成您要查找的内容。
将你的间隔分成更小的间隔。
指定概率密度函数 (PDF)，也可以是一个非常简单的函数，例如阶跃函数。对于高斯分布，您的左右步长将低于中间步长，即（参见下图具有更一般的分布）。
在整个区间内生成一个随机数。如果生成的数字大于您的 PDF 的值，则拒绝生成的数字。
重复这些步骤，直到获得所需的点数
编辑 1
高斯 PDF 上的概念证明。
好的，基本思路如图（a）所示。
- 定义/选择您的概率密度函数 (PDF)。从统计学上讲，PDF 是一个随机变量的函数，它描述了在测量/实验中找到值 x 的概率。一个函数可以是一个随机变量
x 的 PDF，如果它满足：1) f(x) >= 0 和 2) 它被归一化（意味着它求和或积分，直到值 1）。
- 获取 PDF 的最大值 (
max) 和“零点” (z1 < z2)。一些 PDF 的零点可以无穷大。在这种情况下，确定截止点 (z1, z2) 哪个 PDF(z1>x>z2) < eta 您自己选择 eta。基本上意味着，设置一些小的值eta，然后说你的零点是那些PDF(x) 的值小于eta 的值。
- 定义随机发生器的间隔
Ch(z1, z2, max)。这是您生成随机变量的时间间隔。
- 生成一个随机变量 x 使得
z1<x<z2。
- 在
(0, max) 范围内生成第二个不相关的随机变量y。如果y 的值小于PDF(x) 拒绝两个随机生成的值(x,y) 并返回步骤4。如果生成的值y 大于PDF(x) 接受值x 作为随机在分布上生成点并 return 它。
这是重现高斯 PDF 类似行为的代码。
#include "Random.h"
#include <fstream>
using namespace std;
double gaus(double a, double b, double c, double x)
{
return a*exp( -((x-b)*(x-b)/(2*c*c) ));
}
double* random_on_a_gaus_distribution(double inter_a, double inter_b)
{
double res [2];
double a = 1.0; //currently parameters for the Gaussian
double b = 2.0; //are defined here to avoid having
double c = 3.0; //a long function declaration line.
double x = kiss::Ran(inter_a, inter_b);
double y = kiss::Ran(0.0, 1.0);
while (y>gaus(a,b,c,x)) //keep creating values until step 5. is satisfied.
{
x = kiss::Ran(inter_a, inter_b); //this is interval (z1, z2)
y = kiss::Ran(0.0, 1.0); //this is the interval (0, max)
}
res[0] = x;
res[1] = y;
return res; //I return (x,y) for plot reasons, only x is the randomly
} //generated value you're looking for.
void main()
{
double* x;
ofstream f;
f.open("test.txt");
for(int i=0; i<100000; i++)
{
//see bellow how I got -5 and 10 to be my interval (z1, z2)
x = random_on_a_gaus_distribution(-5.0, 10.0);
f << x[0]<<","<<x[1]<<endl;
}
f.close();
}
第 1 步
首先我们在一个名为gaus 的函数中定义高斯PDF 的一般外观。简单的。
然后我们定义一个函数random_on_a_gaus_distribution，它使用一个定义良好的高斯函数。在实验\测量中，我们将通过拟合我们的函数得到系数a, b, c。我为这个例子选择了一些随机的（1、2、3），您可以选择满足您的硬件分配的那些（即：使平均为 7000 的高斯的系数）。
第 2 步和第 3 步
我使用 wolfram mathematica 绘制高斯。使用参数 1,2,3 也可以查看 max 和 (z1, z2) 最合适的值。你可以see the graph yourself。该函数的最大值为 1.0，通过称为 eyeballin 的古老科学方法，我估计截止点为 -5.0 和 10.0。
要使random_on_a_gaus_distribution 更通用，您可以更严格地执行步骤 2）并定义eta，然后在连续点中计算您的函数，直到 PDF 小于 eta。这样做的危险是您的截止点可能相距很远，这对于非常单调的功能可能需要很长时间。此外，您必须自己找到最大值。这通常很棘手，但是一个更简单的问题是最小化函数的负数。对于一般情况，这也可能很棘手，但不是“不可撤销的”。最简单的方法是像我一样作弊，只为几个函数硬编码。
第 4 步和第 5 步
然后你猛扑过去。只要不断创造新的和新的点，直到你达到满意的打击。 请注意返回的数字x是一个随机数。您将无法找到两个连续创建的x 值之间的逻辑链接，或者第一次创建的x 和百万分之一。
但是，在我们分布的 x_max 周围的区间中接受的 x 值的数量大于在 PDF(x) < PDF(x_max) 的区间中创建的 x 值的数量。
这只是意味着您的随机数将在所选区间内以这样的方式加权xi 的值，其中 PDF(xi)<PDF(x)。
我返回了 x 和 y 以便能够绘制下面的图表，但是您想要返回的实际上只是 x。我用 matplotlib 做了图。
最好只显示分布上随机创建的变量的直方图。这表明位于 PDF 函数平均值附近的 x 值最有可能被接受，因此将创建更多具有这些近似值的随机变量。
此外，我假设您会对 Kiss 随机数生成器的实现感兴趣。 你有一个非常好的发电机是非常重要的。我敢说，在某种程度上，吻可能不会削减它（经常使用mersene twister）。
随机.h
#pragma once
#include <stdlib.h>
const unsigned RNG_MAX=4294967295;
namespace kiss{
// unsigned int kiss_z, kiss_w, kiss_jsr, kiss_jcong;
unsigned int RanUns();
void RunGen();
double Ran0(int upper_border);
double Ran(double bottom_border, double upper_border);
}
namespace Crand{
double Ran0(int upper_border);
double Ran(double bottom_border, double upper_border);
}
Kiss.cpp
#include "Random.h"
unsigned int kiss_z = 123456789; //od 1 do milijardu
unsigned int kiss_w = 378295763; //od 1 do milijardu
unsigned int kiss_jsr = 294827495; //od 1 do RNG_MAX
unsigned int kiss_jcong = 495749385; //od 0 do RNG_MAX
//KISS99*
//Autor: George Marsaglia
unsigned int kiss::RanUns()
{
kiss_z=36969*(kiss_z&65535)+(kiss_z>>16);
kiss_w=18000*(kiss_w&65535)+(kiss_w>>16);
kiss_jsr^=(kiss_jsr<<13);
kiss_jsr^=(kiss_jsr>>17);
kiss_jsr^=(kiss_jsr<<5);
kiss_jcong=69069*kiss_jcong+1234567;
return (((kiss_z<<16)+kiss_w)^kiss_jcong)+kiss_jsr;
}
void kiss::RunGen()
{
for (int i=0; i<2000; i++)
kiss::RanUns();
}
double kiss::Ran0(int upper_border)
{
unsigned velicinaIntervala = RNG_MAX / upper_border;
unsigned granicaIzbora= velicinaIntervala*upper_border;
unsigned slucajniBroj = kiss::RanUns();
while(slucajniBroj>=granicaIzbora)
slucajniBroj = kiss::RanUns();
return slucajniBroj/velicinaIntervala;
}
double kiss::Ran (double bottom_border, double upper_border)
{
return bottom_border+(upper_border-bottom_border)*kiss::Ran0(100000)/(100001.0);
}
此外还有标准的 C 随机生成器：
CRands.cpp
#include "Random.h"
//standardni pseudo random generatori iz C-a
double Crand::Ran0(int upper_border)
{
return rand()%upper_border;
}
double Crand::Ran (double bottom_border, double upper_border)
{
return (upper_border-bottom_border)*rand()/((double)RAND_MAX+1);
}
上面的 (b) 图表也值得一提。当您的 PDF 表现非常糟糕时，PDF(x) 在大数字和非常小的数字之间会有很大差异。
问题在于区间区域Ch(x) 将很好地匹配PDF 的极值，但是由于我们也为PDF(x) 的小值创建了一个随机变量y；接受该值的机会微乎其微！此时生成的y 值更有可能总是大于PDF(x)。这意味着您将花费大量周期来创建不会被选中的数字，并且您选择的所有随机数都将在本地绑定到 PDF 的 max。
这就是为什么不要在任何地方都使用相同的Ch(x) 间隔，而是定义一组参数化的间隔通常很有用。然而，这给代码增加了相当多的复杂性。
您在哪里设置限制？边缘案件如何处理？何时以及如何确定您确实需要突然使用这种方法？现在计算max 可能并不那么简单，这取决于您最初设想的方法。
此外，现在您必须纠正这样一个事实，即在您的 Ch(x) 框高度较低的区域更容易接受更多数字，这会使原始 PDF 出现偏差。
这可以通过根据上下边界的高度比对在降低的边界中创建的数字进行加权来纠正，基本上你再重复一次y 步骤。创建一个从 0 到 1 的随机数 z 并将其与比率 lower_height/higher_height 进行比较，保证小于 1。如果z 小于比率：接受x，如果大于则拒绝。
也可以通过编写一个接收对象指针的函数来概括所呈现的代码。通过定义你自己的类，即function，它通常描述函数，在某个点有一个 eval 方法，能够存储你的参数，计算和存储它自己的最大/最小值和零/截止点，你不会有像我一样在函数中传递或定义它们。
祝你好运玩得开心！