??1、隨機數生成：?產生偽隨機數最常用的方法是線性同余法。由線性同余法產生的均勻隨機序列a[1],a[2],...,a[n],...滿足? a[0]=d, a[n]=(b*a[n-1]+c) mod m, n=1,2,...，其中b>=0,c>=0,d>=m。d稱為該隨機序列的種子。如何選取常數b,c和m直接關系到所產生的隨機序列的隨機性能。這是隨機性理論重點研究的內容。從直觀上看，m應該取得充分大。另外應取gcm(m,b)=1，因此可取b為一素數。
? 我們建立一個隨機數類RandomNumber來產生0~(n-1)范圍內的一個隨機數，為了簡化實現，假設我們需要的隨機數不超過16位，即函數的輸入參數n<=65536。在實現時，m取最大數65536，b取一個充分大的素數1194211693，c取12345。RandomNumber類包含一個需由用戶初始化的種子randSeed，給定初始種子后，即可產生與之相應的隨機數序列。randSeed可由用戶指定也可用系統時間自動產生。

[cpp]?view plaincopy

//RandomNumber.hpp:?隨機數類??

#include?<ctime>??

const?unsigned?long?maxshort=65536L;??

const?unsigned?long?multiplier=1194211693L;??

const?unsigned?long?adder=12345L;??

class?RandomNumber{??

private:??

????unsigned?long?randSeed;??//當前種子??

public:??

????RandomNumber(unsigned?long?s=0);??//缺省值0表示由系統自動產生種子??

????unsigned?short?random(unsigned?long?n);?//產生0:n-1之間的隨機整數??

????double?floatRandom(void);??//產生[0,1)之間的隨機實數??

};??

RandomNumber::RandomNumber(unsigned?long?s){??

????if(s==0)??

????????randSeed=time(0);?//用系統時間產生種子??

????else??

????????randSeed=s;?//由用戶提供種子??

}??

unsigned?short?RandomNumber::random(unsigned?long?n){??

????randSeed=multiplier*randSeed+adder;?//用線性同余計算新的種子，因為都是無符號整數，??

????????????????????????????????????????//當結果大于類型上界值時系統會自動進行同余求模運算??

????return?(unsigned?short)((randSeed>>16)%n);??//高16位的隨機性較好，把它映射到0~(n-1)范圍內??

}??

double?RandomNumber::floatRandom(void){??

????return?random(maxshort)/double(maxshort);??

}??

? 函數random(n)在每次計算時，用線性同余式計算出新的種子randSeed，作為產生下一個隨機數的種子。這里的同余運算是由系統自動進行的。因為對于無符號整數的運算，當結果超過unsigned long類型的最大值時，系統會自動進行同余求模運算。新的32位的randSeed是一個隨機數，它的高16位的隨機性比較好，取高16位并映射到0~(n-1)范圍內，就產生了一個我們需要的隨機數。通常random(n)產生的隨機序列是均勻的。函數fRandom()產生[0,1)之間的隨機實數。
??2、用隨機模擬法計算圓周率PI的值：?隨機模擬法的理論基礎是大數定理。設X1,X2,...是相互獨立，服務同一分布的隨機變量序列，且具有數據期望E(Xk)=μ (k,1,2,...)。作前n個變量的算術平均值(X1+X2+...+Xn)/n，則對任意的ε>0，當n趨于無窮大時，有P{ |(X1+X2+...+Xn)/n-μ| < ε }=1。可見，當我們對一個給定的分布進行抽樣，獲得n個樣本X1,X2,...Xn時，樣本的平均值會隨著n的增大依概率收斂于該分布的期望值。這說明我們可以用樣本的均值去估計該分布的期望值，這種隨機模擬方法通常也稱為蒙特卡羅模擬。當樣本個數越多時，估計就越精確。
? 為計算PI的值，用一個邊長為1的單位正方形，里面有一個半徑為1的四分之一圓。向該正方形隨機地投擲（即抽樣）點(x,y)，由于的所投擲的點在正方形內均勻分布，因此它落入圓內的概率為四分之一圓面積與單位正方形面積之比，即PI/4。定義一個隨機變量h(x,y)，當抽樣點(x,y)落入四分之一圓內時其值為1，否則為0，因此h(x,y)的期望值為PI/4。我們可以向正方形內隨機地投擲n個點，以獲得n個樣本h(x1,y1),h(x2,y2),...,h(xn,yn)，設落入四分之圓內的點數為k，則樣本的均值為k/n，用k/n作為PI/4的估計，可計算出PI的近似值，n越大估計就越精確。實現如下：

[cpp]?view plaincopy

#include?"RandomNumber.hpp"??

#include?<iostream>??

//用隨機投點法計算圓周率，n為投的點數??

double?computePI(int?n){??

????static?RandomNumber?dart;??

????int?k=0;??

????for(int?i=1;i<=n;i++){??

????????double?x=dart.floatRandom();??

????????double?y=dart.floatRandom();??

????????if((x*x+y*y)<=1)??

????????????k++;??

????}??

????return?4*k/double(n);??

}??

int?main(void){??

????for(int?i=0;i<10;i++)??//模擬10次??

????????std::cout<<computePI(1000000)<<std::endl;??

????return?0;??

}??

? 模擬10次的運行結果：

[cpp]?view plaincopy

3.14129??

3.14218??

3.14012??

3.14058??

3.1417??

3.14157??

3.14214??

3.13868??

3.14076??

3.13913??

??3、計算定積分：?求函數g(x)在[a,b]上的積分值，記為I。任取一組相互獨立、同分布的隨機變量{Xi}，Xi在[a,b]上服務分布律f(x)（即概率密度函數），其中在{a,b]上f(x)!=0，否則f(x)=0。令h(x)=g(x)/f(x)，其中x<a或x>b時h(x)=0。則{h(Xi)}也是一組相互獨立、同分布的隨機變量。由概率論知，期望值E(h(Xi))等于h(x)f(x)在(-∞,+-∞)上的積分值，即g(x)在[a,b]上的積分值，這正是我們要求的積分值I。因此當我們在[a,b]上抽取分布h(x)的n個樣本h(x1),h(x2),...,h(xn)時，樣本的均值可作為I的估計。我們需要選擇一個有簡便方法可以進行抽樣的概率密度函數f(x)，以得到h(x)。最簡單的構造是f(x)在[a,b]上為一個常數。可取f(x)為[a,b]上的均勻分布，即當x∈[a,b]時f(x)=1/(b-a)，否則f(x)=0。這樣h(x)描述為：當x∈[a,b]時h(x)=(b-a)g(x)，否則h(x)=0。
? 在[a,b]上隨機地投擲n個點，得到n個樣本h(x1),h(x2),...,h(xn)，即(b-a)g(x1),(b-a)g(x2),...,(b-a)g(xn)，用這n個樣本的均值(b-a)[g(x1)+g(x2)+...+g(xn)]/n作為I的估計。實現如下：

[cpp]?view plaincopy

#include?"RandomNumber.hpp"??

#include?<iostream>??

//求函數g(x)在[a,b]上的積分值??

double?integration(double?(*g)(double),double?a,double?b,int?n){??

????static?RandomNumber?rnd;??

????double?x,y=0;??

????for(int?i=1;i<=n;i++){??

????????x=(b-a)*rnd.floatRandom()+a;?//產生[a,b]上的一個隨機點x??

????????y+=g(x);??

????}??

????return?(b-a)*y/(double)n;??

}??

double?func(double?x){??

????return?x;??

}??

int?main(void){??

????for(int?i=0;i<10;i++)??//模擬10次??

????????std::cout<<integration(func,2,4,10000)<<std::endl;??

????return?0;??

}??

? 對g(x)=x，在[2,4]上積分值為6，下面是算法模擬10次的運行結果：

[cpp]?view plaincopy

6.00455??

5.98318??

5.99375??

5.9861??

6.00957??

6.00306??

6.00256??

6.01365??

6.01658??

5.98637??

??4、解非線性方程組：?對一般的非線性方程組fi(x)=0,i=1,2,...,n，其中x為向量x={x1,x2,...,xn}，有許多種數值方法來求解。最常用的有線性化方法和函數極小值方法。我們使用后者，構造目標函數M(x)=[f1(x)]**2+[f2(x)]**2+...+[fn(x)]**2，其中x={x1,x2,...,xn}。即M(x)為各方程函數的平方和。由最優化理論可知，M(x)的極小值點即為所求非線性方程組的一個解。
? 求函數M(x)的極小值點可采用隨機模擬的方法。在指定求根區域D內，選定一個隨機點x0作為隨機搜索的出發點，計算出初始的目標函數值M(x0)。在搜索過程中，假設第j步隨機搜索得到的點為xj，對第j+1步，先計算下一步的隨機搜索方向向量r，然后計算搜索步長a。根據向量r和步長a計算出搜索的步進增量向量dx，由此可得到第j+1步的隨機搜索點xj+1=xj+dx，更新目標函數的值M(xj+1)。如果M(xj+1)小于當前的最小值，表示搜索成功，增加步長a，繼續向前搜索，一旦目標函數值小于給定的充分小的值，搜索結束，返回搜索成功標志，并且最后這一步得到的搜索點為方程組的近似解。如果不小于當前的最小值，則表示當前這次搜索失敗。因為搜索并不總會成功，搜索失敗時得到的解并不是我們需要的解，因此我們必須給定搜索次數的上限，達到上限值搜索結束并返回。根據返回標志是搜索成功還是失敗，來判斷得到的解是否是可行的近似解。實現如下：

[cpp]?view plaincopy

#include?"RandomNumber.hpp"??

#include?<iostream>??

//求非線性方程組的根??

bool?nonlinearEquations(double?(*func[])(double[],int),double*?x0,double*?dx0,double*?x,??

????????double?a0,double?epsilon,double?k,int?n,int?steps,int?M){??

????static?RandomNumber?rnd;??

????bool?success;??//搜索成功標識??

????double?*dx,*r;??

????dx=new?double[n];??//搜索增量向量??

????r=new?double[n];??//搜索方向向量??

????int?mm=0;??//當前搜索失敗次數??

????int?j=0;??//迭代次數??

????double?a=a0;??//步長因子??

????double?fx=0;??//目標函數值??

????for(int?i=0;i<n;i++){??

????????x[i]=x0[i];??//根的初始值??

????????dx[i]=dx0[i];??//搜索增量的初始值??

????}??

????for(int?i=0;i<n;i++)??

????????fx+=func[i](x,n)*func[i](x,n);?//計算初始的目標函數值??

????double?min=fx;??//當前最優值??

????while((min>epsilon)?&&?(j<steps)){??

????????//計算隨機搜索步長??

????????if(fx<min){?//搜索成功??

????????????min=fx;??

????????????a*=k;??//增加步長??

????????????success=true;??

????????}else{??//搜索失敗??

????????????mm++;??

????????????if(mm>M)??a/=k;?//回到原來的步長??

????????????success=false;??

????????}??

????????//計算隨機搜索方向和增量??

????????for(int?i=0;i<n;i++)??

????????????r[i]=2.0*rnd.floatRandom()-1;??//產生搜索方向向量??

????????if(success)??

????????????for(int?i=0;i<n;i++)??

????????????????dx[i]=a*r[i];?//計算搜索增量??

????????else??

????????????for(int?i=0;i<n;i++)??

????????????????dx[i]=a*r[i]-dx[i];??

????????for(int?i=0;i<n;i++)??

????????????x[i]+=dx[i];??//更新方程組的根??

????????//更新目標函數值??

????????fx=0;??

????????for(int?i=0;i<n;i++)??

????????????fx+=func[i](x,n)*func[i](x,n);??

????????j++;??

????}??

????delete[]?dx;??

????delete[]?r;??

????if(fx<=epsilon)??

????????return?true;??

????else??

????????return?false;??

}??

double?func0(double?x[],int?n){??

????return?2*x[0]+x[1]-11;??

}??

double?func1(double?x[],int?n){??

????return?x[0]+2*x[1]-7;??

}??

int?main(void){??

????double?(*func[2])(double?x[],int)={func0,func1};??

????for(int?i=0;i<1000;i++){??//模擬求解1000次??

????????double?x0[2]={5.4,0.8};?//初始解??

????????double?dx0[2]={0.1,0.1};??//初始的步進增量方向??

????????double?x[2]={0,0};??//存放返回的解??

????????if(nonlinearEquations(func,x0,dx0,x,1,0.1,1.5,2,100000,500))??

????????????std::cout<<"x0="<<x[0]<<",?x1="<<x[1]<<std::endl;??

????}??

????return?0;??

}??

? 對簡單的方程組2x+y-11=0,x+2y-7=0，其真實解為(5,1)。算法模擬求解1000次，每次求解的搜索次數為100000，目標函數值的需要達到0.1這么小。運行結果如下：

[cpp]?view plaincopy

x0=5.17093,?x1=0.891772??

x0=5.02852,?x1=0.849973??

x0=5.04405,?x1=0.828128??

x0=5.0859,?x1=1.01767??

x0=5.11071,?x1=0.964664??

?可見求解1000次，只有5次是搜索成功的，得到了可行的近似解。因為算法給定的初始解可能離真實解比較遠，而且搜索方向r是隨機生成的，可能是逼近真實解的方向，也可能是遠離真實解的方向，一旦為后者，會導致目標函數值fx比之前更大，搜索就失敗。另外，算法對給定目標函數值epsilon是很敏感的，如果給定的值比較小，要使每次搜索時的目標函數值都向這個給定值靠近是非常難的，一旦遠離這個值搜索就失敗。可見，搜索失敗的可能性是非常大的。增大給定的目標函數值epsilon可以提高搜索成功的次數，但這樣求得的近似解其精度會降低。
? 在隨機模擬法中，算法實現的關鍵是要進行抽樣，不同的隨機變量分布有不同的抽樣方法。上面介紹的是平面區域或直線區間上滿足均勻分布的點的抽樣，利用[0,1]上均勻分布的隨機數來產生所需的隨機點（即樣本）。一般對均勻分布、指數分布、韋伯分布等，都可以通過直接變換，并利用均勻分布的隨機數來產生所需的樣本，而正態分布、很多離散概率分布（如二項分布、泊松分布）等則需要使用其他的抽樣方法。
??隨機數生成算法：線性同余法
? 隨機模擬法的基本思想：構造一個隨機變量分布使其期望等于問題所求的值、對這個概率分布進行抽樣、用樣本的均值作為問題所求值的估計。

總結

以上是生活随笔為你收集整理的数值概率算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。