隨機數生成方法

時間 2019-11-24

標籤隨機數生成方法简体版

原文原文鏈接

轉自：https://www.cnblogs.com/ECJTUACM-873284962/p/6926203.htmlhtml

一、蒙特卡洛方法

蒙特卡羅方法又稱統計模擬法、隨機抽樣技術，是一種隨機模擬方法，以機率和統計理論方法爲基礎的一種計算方法，是使用隨機數（或更常見的僞隨機數）來解決不少計算問題的方法。將所求解的問題同必定的機率模型相聯繫，用電子計算機實現統計模擬或抽樣，以得到問題的近似解。爲象徵性地代表這一方法的機率統計特徵，數學家馮·諾依曼用聞名世界的賭城——蒙特卡羅命名（就是那個馮·諾依曼）。
蒙特卡羅方法解題過程的主要步驟：
a.針對實際問題創建一個簡單且便於實現的機率統計模型，使所求的量剛好是該模型的機率分佈或數字特徵。
b.對模型的隨機變量創建抽樣方法，在計算機上進行模擬測試，抽取足夠多的隨機數。
c.對模擬實驗結果進行統計分析，給出所求解的「估計」。
d.必要時，改進模型以提升估計精度和減小實驗費用，提升模擬效率。ios

二、馮·諾依曼算法

馮·諾依曼（John von Neumann，1903~1957），20世紀最重要的數學家之一，在現代計算機、博弈論和核武器等諸多領域內有傑出建樹的最偉大的科學全才之一，被稱爲「計算機之父」和「博弈論之父」。主要貢獻是：2進制思想與程序內存思想，固然還有蒙特卡洛方法。經過第一部分，可知，蒙特卡洛方法更多的是一種思想的體現（這點遠不一樣於快排等「嚴格」類算法），下面介紹最多見的一種應用——隨機數生成。dom

三、U(0,1)隨機數的產生函數

對隨機系統進行模擬，便須要產生服從某種分佈的一系列隨機數。最經常使用、最基礎的隨機數是在（0,1）區間內均勻分佈的隨機數，最經常使用的兩類數值計算方法是：乘同餘法和混合同餘法。post

乘同餘法：其中，被稱爲種子，是模，是（0,1）區間的隨機數。測試

混合同餘法：其中，是非負整數。ui

這些隨機數是具備週期性的，模擬參數的選擇不一樣，產生的隨機數質量也有所差別。更復雜的生成方法還有：spa

四、從U(0,1)到其它機率分佈的隨機數code

離散型隨機數的模擬

設隨機變量X的機率分佈爲：，分佈函數有

設隨機變量U~U(0,1)的均勻分佈，則，代表的機率與隨機變量u落在與之間的機率相同。

例如：離散隨機變量X有分佈律

X	0	1	2
P(x)	0.3	0.3	0.4

U是(0,1)的均勻分佈，則有，這樣獲得的x便具備X的分佈律。

連續型隨機變量的模擬

經常使用的有兩種方法：逆變換法和舍選法。逆變換法
定理：設隨機變量Y的分佈函數爲F(y)是連續函數，而U是(0,1)上均勻分佈的隨機變量。另，則X和Y具備相同的分佈。

證實：由定義知，X的分佈函數
因此X和Y具備相同的分佈。
這樣計算得，帶入均勻分佈的U，便可獲得服從的隨機數Y。
例如：設X~U(a,b)，則其分佈函數爲

則。因此生成U(0,1)的隨機數U，則即是來自U(a,b)的隨機數。

有些隨機變量的累計分佈函數不存在或者難以求出，即便存在，但計算困難，因而提出了舍選法
要產生服從的隨機數，設x的值域爲[a,b]，抽樣過程以下：

1.已知隨機分佈且x的取值區間也爲[a,b]，並要求，如圖：

2.從中隨機抽樣得，而後由的均勻分佈抽樣得。
3.接受或捨棄取樣值，若是捨棄該值；返回上一步，不然接受。幾何解釋以下：

常數c的選取：c應該儘量地小，由於抽樣效率與c成反比；通常取。這裏的能夠取均勻分佈，這樣由第二步中兩個均勻分佈便能獲得其餘任意分佈的模擬抽樣。

五、正態隨機數的生成

除了上面的反函數法和舍選法，正態隨機數還能夠根據中心極限定理和Box Muller（座標變換法）獲得。

中心極限定理：若是隨機變量序列獨立同分布，而且具備有限的數學指望和方差，則對於一切有

也就是說，當n個獨立同分布的變量和，服從的正態分佈（n足夠大時）。

設n個獨立同分布的隨機變量，它們服從U(0,1)的均勻分佈，那麼漸近服從正態分佈。

Box Muller方法，設（X,Y）是一對相互獨立的服從正態分佈的隨機變量，則有機率密度函數：

令，其中，則有分佈函數：

令，則分佈函數的反函數得：。

若是服從均勻分佈U(0,1)，則可由模擬生成（也爲均勻分佈，可被代替）。令爲，服從均勻分佈U(0,1)。得：

X和Y均服從正態分佈。用Box Muller方法來生成服從正態分佈的隨機數是十分快捷方便的。

下面介紹幾種簡單的隨機數的算法

1 生成隨機數

通常c語言中提供了隨機數生成函數，

其一是僞隨機數--rand：用於返回一個0-32767之間的僞隨機數；

其二是隨機種子函數--srand：用來初始化隨機數發生器的隨機種子

 1 #include <stdio.h>  2 #include <stdlib.h>  3 #include <time.h>  4  5 int main()  6 {  7 int i,j;  8 srand((int)time(0));  9 for (int i = 0; i < 10; i++) 10  { 11 for (int j = 0; j < 10; j++) 12  { 13 printf("%d ",rand()); 14  } 15 printf("\n"); 16  } 17 return 0; 18 }

固然也能夠生成必定範圍內的隨機數

好比生成0——100之間的隨機數

 1 #include <stdio.h>  2 #include <stdlib.h>  3 #include <time.h>  4  5 int main()  6 {  7 int i,j;  8 srand((int)time(0));  9 for (int i = 0; i < 10; i++) 10  { 11 for (int j = 0; j < 10; j++) 12  { 13 printf("%d ",rand()*100/32767); 14  } 15 printf("\n"); 16  } 17 return 0; 18 }

也能夠生成100——200之間的隨機數

 1 #include <stdio.h>  2 #include <stdlib.h>  3 #include <time.h>  4  5 int main()  6 {  7 int i,j;  8 srand((int)time(0));  9 for (int i = 0; i < 10; i++) 10  { 11 for (int j = 0; j < 10; j++) 12  { 13 printf("%d ",rand()/1000+100); 14  } 15 printf("\n"); 16  } 17 return 0; 18 }

使用rand()函數獲取必定範圍內的一個隨機數

若是想要獲取在必定範圍內的數的話，直接作相應的除法取餘便可。

 1 #include<iostream>  2 #include<ctime>  3 using namespace std;  4 int main()  5 {  6 srand(time(0));  7 for(int i=0;i<10;i++)  8  {  9 //產生10之內的整數 10 cout<<rand()%10<<endl; 11  } 12 }

2 生成[0,1]之間均勻分佈的隨機數算法

在這裏採用一種方式生成隨機數

其中i=1,2,3.。。。

而pi就是地推倒的第i個隨機數

根據經驗，通常選取基數base=256.0,通常爲2的整數倍；另外的兩個常數選取a=17.0 和b=139.0

須要注意

（1）這裏的取模運算是針對浮點型數據的，而c語言中的取模運算不能用於浮點數數據的操做，這樣就須要用戶本身編寫取模的程序；

（2）ri是隨着遞推而每次更新的。所以，若是將這個算法編寫出函數，須要考慮參數是傳值仍是傳地址；

遞推更新，因此在這裏要傳地址，不然得不到結果！

 1 #include <stdio.h>  2  3  4 double rand0_1(double *r)  5 {  6 double base=256.0;  7 double a=17.0;  8 double b=139.0;  9 double temp1=a*(*r)+b; 10 //printf("%lf",temp1); 11 double temp2=(int)(temp1/base); //獲得餘數 12 double temp3=temp1-temp2*base; 13 //printf("%lf\n",temp2); 14 //printf("%lf\n",temp3); 15 *r=temp3; 16 double p=*r/base; 17 return p; 18 } 19 20 int main() 21 { 22 double r=5.0; 23 printf("output 10 number between 0 and 1:\n"); 24 for (int i = 0; i < 10; i++) 25  { 26 printf("%10.5lf\n",rand0_1(&r)); 27  } 28 return 0; 29 }

3 產生任意範圍內的隨機數，好比產生[m,n]之間的隨機數

這個很容易，只要將以前的[0,1]之間的隨機數這樣處理就好了

m+(m-n)*rand0_1(&r)就好了；

 1 #include <stdio.h>  2  3  4 double rand0_1(double *r)  5 {  6 double base=256.0;  7 double a=17.0;  8 double b=139.0;  9 double temp1=a*(*r)+b; 10 //printf("%lf",temp1); 11 double temp2=(int)(temp1/base); //獲得餘數 12 double temp3=temp1-temp2*base; 13 //printf("%lf\n",temp2); 14 //printf("%lf\n",temp3); 15 *r=temp3; 16 double p=*r/base; 17 return p; 18 } 19 20 int main() 21 { 22 double m=1.0,n=5.0; 23 double r=5.0; 24 printf("output 10 number between 0 and 1:\n"); 25 for (int i = 0; i < 10; i++) 26  { 27 printf("%10.5lf\n",m+(n-m)*rand0_1(&r)); 28  } 29 return 0; 30 }

4 正態分佈的隨機數生成算法

符合正太分佈的隨機數在研究中也很重要，下面給出一種生成正態分佈數的方法

其中Ri表示[0,1]之間均勻分佈的隨機數；

u爲均值，爲方差，當n趨向於無窮大的時候，獲得隨機的隨機分佈爲正態分佈；

 1 #include <stdio.h>  2 #include <math.h>  3  4 double rand0_1(double *r)  5 {  6 double base=256.0;  7 double a=17.0;  8 double b=139.0;  9 double temp1=a*(*r)+b; 10 //printf("%lf",temp1); 11 double temp2=(int)(temp1/base); //獲得餘數 12 double temp3=temp1-temp2*base; 13 //printf("%lf\n",temp2); 14 //printf("%lf\n",temp3); 15 *r=temp3; 16 double p=*r/base; 17 return p; 18 } 19 20 double random_normality(double u,double t,double *r ,double n) 21 { 22 double total=0.0; 23 double result; 24 for (int i = 0; i < n; i++) 25  { 26 total+=rand0_1(r); 27  } 28 result=u+t*(total-n/2)/sqrt(n/12); 29 return result; 30 } 31 32 int main() 33 { 34 double r=5.0; 35 double u=2.0; 36 double t=3.5; 37 double n=12; 38 printf("output 10 number between 0 and 1:\n"); 39 for (int i = 0; i < 10; i++) 40  { 41 printf("%10.5lf\n",random_normality(u,t,&r,n)); 42  } 43 return 0; 44 }

補充知識點：leveldb中使用了一個簡單的方式來實現隨機化數；算法的核心是seed_ = (seed_ * A) % M,

下面把源代碼貼出來，不難，能夠和上面的參考下

 1 private:  2  uint32_t seed_;  3 public:  4 explicit Random(uint32_t s) : seed_(s & 0x7fffffffu) {  5 // Avoid bad seeds.  6 if (seed_ == 0 || seed_ == 2147483647L) {  7 seed_ = 1;  8  }  9  } 10  uint32_t Next() { 11 static const uint32_t M = 2147483647L; // 2^31-1 12 static const uint64_t A = 16807; // bits 14, 8, 7, 5, 2, 1, 0 13 // We are computing 14 // seed_ = (seed_ * A) % M, where M = 2^31-1 15 // 16 // seed_ must not be zero or M, or else all subsequent computed values 17 // will be zero or M respectively. For all other values, seed_ will end 18 // up cycling through every number in [1,M-1] 19 uint64_t product = seed_ * A; 20 21 // Compute (product % M) using the fact that ((x << 31) % M) == x. 22 seed_ = static_cast<uint32_t>((product >> 31) + (product & M)); 23 // The first reduction may overflow by 1 bit, so we may need to 24 // repeat. mod == M is not possible; using > allows the faster 25 // sign-bit-based test. 26 if (seed_ > M) { 27 seed_ -= M; 28  } 29 return seed_; 30  } 31 // Returns a uniformly distributed value in the range [0..n-1] 32 // REQUIRES: n > 0 33 uint32_t Uniform(int n) { return Next() % n; } 34 35 // Randomly returns true ~"1/n" of the time, and false otherwise. 36 // REQUIRES: n > 0 37 bool OneIn(int n) { return (Next() % n) == 0; } 38 39 // Skewed: pick "base" uniformly from range [0,max_log] and then 40 // return "base" random bits. The effect is to pick a number in the 41 // range [0,2^max_log-1] with exponential bias towards smaller numbers. 42 uint32_t Skewed(int max_log) { 43 return Uniform(1 << Uniform(max_log + 1)); 44  } 45 };