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); //用系统时间产生种子
- 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越大估计就越精确。实现如下：
- #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.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的估计。实现如下：
- #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次的运行结果：
- 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，继续向前搜索，一旦目标函数值小于给定的充分小的值，搜索结束，返回搜索成功标志，并且最后这一步得到的搜索点为方程组的近似解。如果不小于当前的最小值，则表示当前这次搜索失败。因为搜索并不总会成功，搜索失败时得到的解并不是我们需要的解，因此我们必须给定搜索次数的上限，达到上限值搜索结束并返回。根据返回标志是搜索成功还是失败，来判断得到的解是否是可行的近似解。实现如下：
- #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这么小。运行结果如下：
- 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]上均匀分布的随机数来产生所需的随机点（即样本）。一般对均匀分布、指数分布、韦伯分布等，都可以通过直接变换，并利用均匀分布的随机数来产生所需的样本，而正态分布、很多离散概率分布（如二项分布、泊松分布）等则需要使用其他的抽样方法。
随机数生成算法：线性同余法
随机模拟法的基本思想：构造一个随机变量分布使其期望等于问题所求的值、对这个概率分布进行抽样、用样本的均值作为问题所求值的估计。
扫一扫
2923
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
