数值概率算法

  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可由用户指定也可用系统时间自动产生。

//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); //用系统时间产生种子elserandSeed=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越大估计就越精确。实现如下：

#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次的运行结果：

3.141293.142183.140123.140583.14173.141573.142143.138683.140763.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的估计。实现如下：

#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]上的一个随机点xy+=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次的运行结果：

6.004555.983185.993755.98616.009576.003066.002566.013656.016585.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，继续向前搜索，一旦目标函数值小于给定的充分小的值，搜索结束，返回搜索成功标志，并且最后这一步得到的搜索点为方程组的近似解。如果不小于当前的最小值，则表示当前这次搜索失败。因为搜索并不总会成功，搜索失败时得到的解并不是我们需要的解，因此我们必须给定搜索次数的上限，达到上限值搜索结束并返回。根据返回标志是搜索成功还是失败，来判断得到的解是否是可行的近似解。实现如下：

#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]; //计算搜索增量elsefor(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;elsereturn 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这么小。运行结果如下：

x0=5.17093, x1=0.891772x0=5.02852, x1=0.849973x0=5.04405, x1=0.828128x0=5.0859, x1=1.01767x0=5.11071, x1=0.964664

 可见求解1000次，只有5次是搜索成功的，得到了可行的近似解。因为算法给定的初始解可能离真实解比较远，而且搜索方向r是随机生成的，可能是逼近真实解的方向，也可能是远离真实解的方向，一旦为后者，会导致目标函数值fx比之前更大，搜索就失败。另外，算法对给定目标函数值epsilon是很敏感的，如果给定的值比较小，要使每次搜索时的目标函数值都向这个给定值靠近是非常难的，一旦远离这个值搜索就失败。可见，搜索失败的可能性是非常大的。增大给定的目标函数值epsilon可以提高搜索成功的次数，但这样求得的近似解其精度会降低。
  在随机模拟法中，算法实现的关键是要进行抽样，不同的随机变量分布有不同的抽样方法。上面介绍的是平面区域或直线区间上满足均匀分布的点的抽样，利用[0,1]上均匀分布的随机数来产生所需的随机点（即样本）。一般对均匀分布、指数分布、韦伯分布等，都可以通过直接变换，并利用均匀分布的随机数来产生所需的样本，而正态分布、很多离散概率分布（如二项分布、泊松分布）等则需要使用其他的抽样方法。
  随机数生成算法：线性同余法
  随机模拟法的基本思想：构造一个随机变量分布使其期望等于问题所求的值、对这个概率分布进行抽样、用样本的均值作为问题所求值的估计。