数值概率算法

  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 plain copy
  1. //RandomNumber.hpp: 随机数类  
  2. #include <ctime>  
  3. const unsigned long maxshort=65536L;  
  4. const unsigned long multiplier=1194211693L;  
  5. const unsigned long adder=12345L;  
  6. class RandomNumber{  
  7. private:  
  8.     unsigned long randSeed;  //当前种子  
  9. public:  
  10.     RandomNumber(unsigned long s=0);  //缺省值0表示由系统自动产生种子  
  11.     unsigned short random(unsigned long n); //产生0:n-1之间的随机整数  
  12.     double floatRandom(void);  //产生[0,1)之间的随机实数  
  13. };  
  14. RandomNumber::RandomNumber(unsigned long s){  
  15.     if(s==0)  
  16.         randSeed=time(0); //用系统时间产生种子  
  17.     else  
  18.         randSeed=s; //由用户提供种子  
  19. }  
  20. unsigned short RandomNumber::random(unsigned long n){  
  21.     randSeed=multiplier*randSeed+adder; //用线性同余计算新的种子,因为都是无符号整数,  
  22.                                         //当结果大于类型上界值时系统会自动进行同余求模运算  
  23.     return (unsigned short)((randSeed>>16)%n);  //高16位的随机性较好,把它映射到0~(n-1)范围内  
  24. }  
  25. double RandomNumber::floatRandom(void){  
  26.     return random(maxshort)/double(maxshort);  
  27. }  

  函数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 plain copy
  1. #include "RandomNumber.hpp"  
  2. #include <iostream>  
  3. //用随机投点法计算圆周率,n为投的点数  
  4. double computePI(int n){  
  5.     static RandomNumber dart;  
  6.     int k=0;  
  7.     for(int i=1;i<=n;i++){  
  8.         double x=dart.floatRandom();  
  9.         double y=dart.floatRandom();  
  10.         if((x*x+y*y)<=1)  
  11.             k++;  
  12.     }  
  13.     return 4*k/double(n);  
  14. }  
  15. int main(void){  
  16.     for(int i=0;i<10;i++)  //模拟10次  
  17.         std::cout<<computePI(1000000)<<std::endl;  
  18.     return 0;  
  19. }  

  模拟10次的运行结果:

[cpp]  view plain copy
  1. 3.14129  
  2. 3.14218  
  3. 3.14012  
  4. 3.14058  
  5. 3.1417  
  6. 3.14157  
  7. 3.14214  
  8. 3.13868  
  9. 3.14076  
  10. 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 plain copy
  1. #include "RandomNumber.hpp"  
  2. #include <iostream>  
  3. //求函数g(x)在[a,b]上的积分值  
  4. double integration(double (*g)(double),double a,double b,int n){  
  5.     static RandomNumber rnd;  
  6.     double x,y=0;  
  7.     for(int i=1;i<=n;i++){  
  8.         x=(b-a)*rnd.floatRandom()+a; //产生[a,b]上的一个随机点x  
  9.         y+=g(x);  
  10.     }  
  11.     return (b-a)*y/(double)n;  
  12. }  
  13. double func(double x){  
  14.     return x;  
  15. }  
  16. int main(void){  
  17.     for(int i=0;i<10;i++)  //模拟10次  
  18.         std::cout<<integration(func,2,4,10000)<<std::endl;  
  19.     return 0;  
  20. }  

  对g(x)=x,在[2,4]上积分值为6,下面是算法模拟10次的运行结果:

[cpp]  view plain copy
  1. 6.00455  
  2. 5.98318  
  3. 5.99375  
  4. 5.9861  
  5. 6.00957  
  6. 6.00306  
  7. 6.00256  
  8. 6.01365  
  9. 6.01658  
  10. 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 plain copy
  1. #include "RandomNumber.hpp"  
  2. #include <iostream>  
  3. //求非线性方程组的根  
  4. bool nonlinearEquations(double (*func[])(double[],int),double* x0,double* dx0,double* x,  
  5.         double a0,double epsilon,double k,int n,int steps,int M){  
  6.     static RandomNumber rnd;  
  7.     bool success;  //搜索成功标识  
  8.     double *dx,*r;  
  9.     dx=new double[n];  //搜索增量向量  
  10.     r=new double[n];  //搜索方向向量  
  11.     int mm=0;  //当前搜索失败次数  
  12.     int j=0;  //迭代次数  
  13.     double a=a0;  //步长因子  
  14.     double fx=0;  //目标函数值  
  15.     for(int i=0;i<n;i++){  
  16.         x[i]=x0[i];  //根的初始值  
  17.         dx[i]=dx0[i];  //搜索增量的初始值  
  18.     }  
  19.     for(int i=0;i<n;i++)  
  20.         fx+=func[i](x,n)*func[i](x,n); //计算初始的目标函数值  
  21.     double min=fx;  //当前最优值  
  22.     while((min>epsilon) && (j<steps)){  
  23.         //计算随机搜索步长  
  24.         if(fx<min){ //搜索成功  
  25.             min=fx;  
  26.             a*=k;  //增加步长  
  27.             success=true;  
  28.         }else{  //搜索失败  
  29.             mm++;  
  30.             if(mm>M)  a/=k; //回到原来的步长  
  31.             success=false;  
  32.         }  
  33.         //计算随机搜索方向和增量  
  34.         for(int i=0;i<n;i++)  
  35.             r[i]=2.0*rnd.floatRandom()-1;  //产生搜索方向向量  
  36.         if(success)  
  37.             for(int i=0;i<n;i++)  
  38.                 dx[i]=a*r[i]; //计算搜索增量  
  39.         else  
  40.             for(int i=0;i<n;i++)  
  41.                 dx[i]=a*r[i]-dx[i];  
  42.         for(int i=0;i<n;i++)  
  43.             x[i]+=dx[i];  //更新方程组的根  
  44.         //更新目标函数值  
  45.         fx=0;  
  46.         for(int i=0;i<n;i++)  
  47.             fx+=func[i](x,n)*func[i](x,n);  
  48.         j++;  
  49.     }  
  50.     delete[] dx;  
  51.     delete[] r;  
  52.     if(fx<=epsilon)  
  53.         return true;  
  54.     else  
  55.         return false;  
  56. }  
  57. double func0(double x[],int n){  
  58.     return 2*x[0]+x[1]-11;  
  59. }  
  60. double func1(double x[],int n){  
  61.     return x[0]+2*x[1]-7;  
  62. }  
  63. int main(void){  
  64.     double (*func[2])(double x[],int)={func0,func1};  
  65.     for(int i=0;i<1000;i++){  //模拟求解1000次  
  66.         double x0[2]={5.4,0.8}; //初始解  
  67.         double dx0[2]={0.1,0.1};  //初始的步进增量方向  
  68.         double x[2]={0,0};  //存放返回的解  
  69.         if(nonlinearEquations(func,x0,dx0,x,1,0.1,1.5,2,100000,500))  
  70.             std::cout<<"x0="<<x[0]<<", x1="<<x[1]<<std::endl;  
  71.     }  
  72.     return 0;  
  73. }  

  对简单的方程组2x+y-11=0,x+2y-7=0,其真实解为(5,1)。算法模拟求解1000次,每次求解的搜索次数为100000,目标函数值的需要达到0.1这么小。运行结果如下:

[cpp]  view plain copy
  1. x0=5.17093, x1=0.891772  
  2. x0=5.02852, x1=0.849973  
  3. x0=5.04405, x1=0.828128  
  4. x0=5.0859, x1=1.01767  
  5. x0=5.11071, x1=0.964664  

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值