随机数生成算法

时间：2024-11-10 20:22:31

转自：/ECJTUACM-873284962/p/

1、蒙特卡罗法

蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，数学家冯·诺依曼用闻名世界的赌城——蒙特卡罗命名（就是那个冯·诺依曼）。
蒙特卡罗方法解题过程的主要步骤：
a.针对实际问题建立一个简单且便于实现的概率统计模型，使所求的量恰好是该模型的概率分布或数字特征。
b.对模型的随机变量建立抽样方法，在计算机上进行模拟测试，抽取足够多的随机数。
c.对模拟实验结果进行统计分析，给出所求解的“估计”。
d.必要时，改进模型以提高估计精度和减少实验费用，提高模拟效率。

2、冯▪诺依曼

冯·诺依曼（John von Neumann，1903~1957），20世纪最重要的数学家之一，在现代计算机、博弈论和核武器等诸多领域内有杰出建树的最伟大的科学全才之一，被称为“计算机之父”和“博弈论之父”。主要贡献是：二进制思想与程序内存思想，当然还有蒙特卡洛方法。通过第一部分，可知，蒙特卡洛方法更多的是一种思想的体现（这点远不同于快排等“严格”类算法），下面介绍最常见的一种应用——随机数生成。

3、U(0,1)随机数的产生

对随机系统进行模拟，便需要产生服从某种分布的一系列随机数。最常用、最基础的随机数是在（0,1）区间内均匀分布的随机数，最常用的两类数值计算方法是：乘同余法和混合同余法。
乘同余法：
　这里写图片描述
其中， x0 x 0 被称为种子， M M 是模，rn是（0,1）区间的随机数。

混合同余法：
　这里写图片描述
其中，C是非负整数。
这些随机数是具有周期性的，模拟参数的选择不同，产生的随机数质量也有所差异。更复杂的生成方法还有：
　这里写图片描述

4、从U(0,1)到其它概率分布的随机数

（ａ）离散型随机数的模拟
设随机变量X的概率分布为： p(X=xi)=pi p ( X = x i ) = p i ，分布函数有 P(0)=0,P(n)=∑pi P ( 0 ) = 0 , P ( n ) = ∑ p i ，设随机变量U~U(0,1)的均匀分布，则 p(P(n−1)<u<P(n)−P(n−1)=p(n)) p ( P ( n − 1 ) < u < P ( n ) − P ( n − 1 ) = p ( n ) ) 表明 X=xn(即pn) X = x n ( 即 p n ) 的概率与随机变量u落在 P(n−1) P ( n − 1 ) 与 P(n) P ( n ) 之间的概率相同。

例如：离散随机变量X有分布律

X 0 1 2
P(x) 0.3 0.3 0.4

U是(0,1)的均匀分布，则有
　clip_image002[28]
这样得到的x便具有X的分布律。
（ｂ）连续型随机变量的模拟

常用的有两种方法：逆变换法和舍选法。
逆变换法：设随机变量Y的分布函数为F(y)是连续函数，而U是(0,1)上均匀分布的随机变量。另 X=F−1(U) X = F − 1 ( U ) ，则X和Y具有相同的分布。

证明：由定义知，X的分布函数：
clip_image002[32]
所以X和Y具有相同的分布。这样计算得 F−1(∙) F − 1 ( • ) ，带入均匀分布的U，即可得到服从 f(∙) f ( • ) 的随机数Y。

例如：设X~U(a,b)，则其分布函数为：
　clip_image002[42]
则
clip_image002[44]
所以生成U(0,1)的随机数U，则 a+(b−a)u a + ( b − a ) u 便是来自U(a,b)的随机数。有些随机变量的累计分布函数不存在或者难以求出，即使存在，但计算困难，于是提出了
舍选法：要产生服从 p(x) p ( x ) 的随机数，设x的值域为[a,b]，抽样过程如下：
1、已知随机分布 Q(x) Q ( x ) 且x的取值区间也为[a,b]，并要求 c×Q(x)﹥p(x)，xϵ[a,b] c × Q ( x ) ﹥ p ( x ) ， x ϵ [ a , b ] ，如图
这里写图片描述
2、从 Q(x) Q ( x ) 中随机抽样得 x∗ x ∗ ，然后由 U(0,c×Q(x∗)) U ( 0 , c × Q ( x ∗ ) ) 的均匀分布抽样得 u u 。
3、接受或舍弃取样值x∗，如果 u>P(x∗) u > P ( x ∗ ) 舍弃该值；返回上一步，否则接受。几何解释如下：
这里写图片描述
常数c的选取：c应该尽可能地小，因为抽样效率与c成反比；一般取 c=maxP(x)/Q(x),xϵ[a,b] c = m a x P ( x ) / Q ( x ) , x ϵ [ a , b ] 。这里的 Q(x) Q ( x ) 可以取均匀分布，这样由第二种中两个均匀分布便能得到其他任意分布的模拟抽样。

5、正态随机数的生成

除了上面的反函数法和舍选法，正态随机数还可以根据中心极限定理和Box Muller（坐标变换法）得到。
中心极限定理：如果随机变量序列 X1,X2,...,Xn,.. X 1 , X 2 , . . . , X n , . . 独立同分布，并且具有有限的数学期望和方差 E(Zi)=μ,D(Xi)=σ2(i=1,2,...) E ( Z i ) = μ , D ( X i ) = σ 2 ( i = 1 , 2 , . . . ) ，则对于一切x∈R有
clip_image002[80]
也就是说，当n个独立同分布的变量和，服从 N(nμ,nσ2) N ( n μ , n σ 2 ) 的正态分布（n足够大时）。
设n个独立同分布的随机变量 U1,U2,...,Un U 1 , U 2 , . . . , U n ，它们服从U(0,1)的均匀分布，那么
clip_image002[86]
渐近服从正态分布 N(0,1) N ( 0 , 1 ) 。
Box Muller方法，设（X,Y）是一对相互独立的服从正态分布 N(0,1) N ( 0 , 1 ) 的随机变量，则有概率密度函数：
clip_image002[90]
令 x=Rcosθ,y=Rsinθ x = R c o s θ , y = R s i n θ ，其中 x∈(0,2π) x ∈ ( 0 , 2 π ) ，则R有分布函数：
　clip_image002[99]
令 FR(r)=1−e−r22 F R ( r ) = 1 − e − r 2 2 ，则分布函数的反函数得：
　clip_image002[103]
如果 U1 U 1 服从均匀分布U(0,1)，则R可由 −2lnU1−−−−−−√ − 2 l n U 1 模拟生成（ (1−U1) ( 1 − U 1 ) 也为均匀分布，可被 U1 U 1 代替）。令 θ θ 为 2πU2 2 π U 2 ， U2 U 2 服从均匀分布U(0,1)。得：
　clip_image002[124]
X和Y均服从正态分布。用Box Muller方法来生成服从正态分布的随机数是十分快捷方便的。

下面介绍几种简单的随机数的算法

1、 生成随机数
一般c语言中提供了随机数生成函数，
其一是伪随机数–rand：用于返回一个0-32767之间的伪随机数；
其二是随机种子函数–srand：用来初始化随机数发生器的随机种子

#include <>
#include <>
#include <>

int main()
{
    int i,j;
    srand((int)time(0));
    for (int i = 0; i < 10; i++)
    {
        for (int j = 0; j < 10; j++)
        {
            printf("%d  ",rand());
        }
        printf("\n");
    }
    return 0;
}

当然也可以生成一定范围内的随机数，比如生成0——100之间的随机数。

#include <>
#include <>
#include <>

int main()
{
    int i,j;
    srand((int)time(0));
    for (int i = 0; i < 10; i++)
    {
        for (int j = 0; j < 10; j++)
        {
            printf("%d  ",rand()*100/32767);
        }
        printf("\n");
    }
    return 0;
}

也可以生成100——200之间的随机数

#include <>
#include <>
#include <>

int main()
{
    int i,j;
    srand((int)time(0));
    for (int i = 0; i < 10; i++)
    {
        for (int j = 0; j < 10; j++)
        {
            printf("%d  ",rand()/1000+100);
        }
        printf("\n");
    }
    return 0;
}

使用rand()函数获取一定范围内的一个随机数，如果想要获取在一定范围内的数的话，直接做相应的除法取余即可。

#include<iostream>
#include<ctime>
using namespace std;
int main()
{
 srand(time(0));
 for(int i=0;i<10;i++)
 {
  //产生10以内的整数
  cout<<rand()%10<<endl;
 }
}

2 、生成[0,1]之间均匀分布的随机数算法

ri=mod(a∗ri−1+b,base) r i = m o d ( a ∗ r i − 1 + b , b a s e )
pi=ri/base p i = r i / b a s e

在这里采用第二种方式生成随机数。其中i=1,2,3.。。。而pi就是地推倒的第i个随机数，根据经验，一般选取基数base=256.0,一般为2的整数倍；另外的两个常数选取a=17.0 和b=139.0，需要注意：
（1）这里的取模运算是针对浮点型数据的，而c语言中的取模运算不能用于浮点数数据的操作，这样就需要用户自己编写取模的程序；
（2）ri是随着递推而每次更新的。因此，如果将这个算法编写出函数，需要考虑参数是传值还是传地址；
递推更新，所以在这里要传地址，否则得不到结果！

#include <>


double rand0_1(double *r)
{
    double base=256.0;
    double a=17.0;
    double b=139.0;
    double temp1=a*(*r)+b;
    //printf("%lf",temp1);
    double temp2=(int)(temp1/base); //得到余数
    double temp3=temp1-temp2*base;
    //printf("%lf\n",temp2);
    //printf("%lf\n",temp3);
    *r=temp3;
    double p=*r/base;
    return p;
}

int main()
{
    double r=5.0;
    printf("output 10 number between 0 and 1:\n");
    for (int i = 0; i < 10; i++)
    {
        printf("%10.5lf\n",rand0_1(&r));
    }
    return 0;
}

3 、产生任意范围内的随机数
比如产生[m,n]之间的随机数，这个很容易，只要将之前的[0,1]之间的随机数这样处理就行了，m+(m-n)*rand0_1(&r)就行了。

#include <>


double rand0_1(double *r)
{
    double base=256.0;
    double a=17.0;
    double b=139.0;
    double temp1=a*(*r)+b;
    //printf("%lf",temp1);
    double temp2=(int)(temp1/base); //得到余数
    double temp3=temp1-temp2*base;
    //printf("%lf\n",temp2);
    //printf("%lf\n",temp3);
    *r=temp3;
    double p=*r/base;
    return p;
}

int main()
{
    double m=1.0,n=5.0;
    double r=5.0;
    printf("output 10 number between 0 and 1:\n");
    for (int i = 0; i < 10; i++)
    {
        printf("%10.5lf\n",m+(n-m)*rand0_1(&r));
    }
    return 0;
}

4、正态分布的随机数生成算法
符合正太分布的随机数在研究中也很重要，下面给出一种生成正态分布数的方法

random_normality=u+σ(∑Ri)−n/2n/12−−−−√ r a n d o m _ n o r m a l i t y = u + σ ( ∑ R i ) − n / 2 n / 12

其中Ri表示[0,1]之间均匀分布的随机数；u为均值， σ2 σ 2 为方差，当n趋向于无穷大的时候，得到随机的随机分布为正态分布；

#include <>
#include <>

double rand0_1(double *r)
{
      double base=256.0;
      double a=17.0;
      double b=139.0;
      double temp1=a*(*r)+b;
      //printf("%lf",temp1);
      double temp2=(int)(temp1/base); //得到余数
      double temp3=temp1-temp2*base;
      //printf("%lf\n",temp2);
      //printf("%lf\n",temp3);
      *r=temp3;
      double p=*r/base;
      return p;
}

double random_normality(double u,double t,double *r ,double n)
{
      double total=0.0;
      double result;
      for (int i = 0; i < n; i++)
      {
            total+=rand0_1(r);
      }
      result=u+t*(total-n/2)/sqrt(n/12);
      return result;
}

int main()
{
      double r=5.0;
      double u=2.0;
      double t=3.5;
      double n=12;
      printf("output 10 number between 0 and 1:\n");
      for (int i = 0; i < 10; i++)
      {
            printf("%10.5lf\n",random_normality(u,t,&r,n));
      }
      return 0;
}

补充知识点：leveldb中使用了一个简单的方式来实现随机化数；算法的核心是seed_ =(seed_ * A) % M，下面把源代码贴出来，不难，可以和上面的参考下

private:
  uint32_t seed_;
 public:
  explicit Random(uint32_t s) : seed_(s & 0x7fffffffu) {
    // Avoid bad seeds.
    if (seed_ == 0 || seed_ == 2147483647L) {
      seed_ = 1;
    }
  }
  uint32_t Next() {
    static const uint32_t M = 2147483647L;   // 2^31-1
    static const uint64_t A = 16807;  // bits 14, 8, 7, 5, 2, 1, 0
    // We are computing
    //       seed_ = (seed_ * A) % M,    where M = 2^31-1
    //
    // seed_ must not be zero or M, or else all subsequent computed values
    // will be zero or M respectively.  For all other values, seed_ will end
    // up cycling through every number in [1,M-1]
    uint64_t product = seed_ * A;

    // Compute (product % M) using the fact that ((x << 31) % M) == x.
    seed_ = static_cast<uint32_t>((product >> 31) + (product & M));
    // The first reduction may overflow by 1 bit, so we may need to
    // repeat.  mod == M is not possible; using > allows the faster
    // sign-bit-based test.
    if (seed_ > M) {
      seed_ -= M;
    }
    return seed_;
  }
  // Returns a uniformly distributed value in the range [0..n-1]
  // REQUIRES: n > 0
  uint32_t Uniform(int n) { return Next() % n; }

  // Randomly returns true ~"1/n" of the time, and false otherwise.
  // REQUIRES: n > 0
  bool OneIn(int n) { return (Next() % n) == 0; }

  // Skewed: pick "base" uniformly from range [0,max_log] and then
  // return "base" random bits.  The effect is to pick a number in the
  // range [0,2^max_log-1] with exponential bias towards smaller numbers.
  uint32_t Skewed(int max_log) {
    return Uniform(1 << Uniform(max_log + 1));
  }
};

总之，做个简单的总结
C语言/C++怎样产生随机数：这里要用到的是rand()函数, srand()函数,和time()函数。
需要说明的是，iostream头文件中就有srand函数的定义，不需要再额外引入;而使用time()函数需要引入ctime头文件。
使用rand()函数获取一个随机数
如果你只要产生随机数而不需要设定范围的话，你只要用rand()就可以了：rand()会返回一随机数值, 范围在0至RAND_MAX 间。RAND_MAX定义在, 其值为2147483647。
使用rand函数和time函数
我们上面已经可以获取随机数了，为什么还需要使用time函数呢？我们通过多次运行发现，该程序虽然生成了10个随机数，但是这个10个随机数是固定的，也就是说并不随着时间的变化而变化。
这与srand()函数有关。srand()用来设置rand()产生随机数时的随机数种子。在调用rand()函数产生随机数前，必须先利用srand()设好随机数种子（seed）, 如果未设随机数种子, rand()在调用时会自动设随机数种子为1。
上面的例子就是因为没有设置随机数种子，每次随机数种子都自动设成相同值1 ，进而导致rand()所产生的随机数值都一样。
srand()函数定义 ： void srand (unsigned int seed);
通常可以利用geypid()或time(0)的返回值来当做seed
如果你用time(0)的话，要加入头文件#include
time(0)或者time(NULL)返回的是系统的时间（从1970.1.1午夜算起），单位：秒

作者： Angel_Kitty
出处：/ECJTUACM-873284962/
版权声明：本文版权归作者和博客园共有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文链接。<

相关文章