求助?关于蒙特卡罗算法
谁有相关个实例,小弟感激不尽 <P>同问，也想学习一下</P> 这个问题我也想研究一个，希望大虾们指教! 例如计算平面坐标系里曲线与坐标所围的面积,用积分可以用<STRONG>蒙特卡罗算法也可以,说白了就是生成很多个点记录下在曲线先的点数和一共生成的点数,两者之比.为了计算方便一般会有一个可比较的面积了</STRONG><BR> <P>这个我也不大清楚，谁有相关资料啊？</P> <P>偶也想研究下~~~~~~~~~</P> <table bordercolor="#aecaae" cellspacing="0" cellpadding="0" width="100%" align="center" border="1"><tbody><tr><td valign="bottom" bordercolor="#ffffff" align="center" width="100%" height="60" style="FONT-SIZE: 18px;"><strong><font face="楷体_GB2312" color="#003399" size="5"><b>集合相等的蒙特卡罗算法<br/></b></font></strong></td></tr><tr><td bordercolor="#ffffff" align="center" width="100%"><hr align="center" width="100%" noshade="true" size="1"/></td></tr><tr><td bordercolor="#ffffff" align="center" width="100%" style="FONT-SIZE: 9pt;">发表日期：2006年1月12日 已经有66位读者读过此文</td></tr><tr><td bordercolor="#ffffff" align="center" width="100%"><table cellspacing="0" cellpadding="0" width="98%" align="center" border="0" style="TABLE-LAYOUT: fixed;"><tbody><tr><td align="center" width="100%"></td></tr><tr><td style="WORD-WRAP: break-word;"><font class="news"><br/><pre><font size="3" style="BACKGROUND-COLOR: #f0f8ff;">#include<iostream.h>#include<fstream.h>
#include<math.h>
#include<time.h>
//============随机数类=================
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 fRandom(void); //产生[0，1)之间的随机实数
};
RandomNumber::RandomNumber(unsigned long s)
{//产生种子
if(s==0)
randSeed=time(0); //用系统时间产生种子
else
randSeed=s; //由用户提供种子
}
unsigned short RandomNumber::Random(unsigned long n)
{//产生0:n-1之间的随机整数
randSeed = multiplier * randSeed + adder;
return(unsigned short)((randSeed>>16) % n);
}
double RandomNumber::fRandom(void)
{//产生[0，1)之间的随机实数
return Random(maxshort)/double(maxshort);
}
//===================================================
//=============合并排序算法====================
template <class T>
void Merge(T *c,T *d,int l,int m,int r)
{
int i=l,
j=m+1,
k=l;
while((i<=m)&&(j<=r))
if (c <=c) d=c;
else d=c;
if(i>m) for (int q=j;q<=r;q++)
d=c;
else for(int q=i;q<=m;q++)
d=c;
}
template <class T>
void MergePass(T *x,T *y,int s,int n)
{
int i=0;
while(i<=n-2*s)
{
Merge(x,y,i,i+s-1,i+2*s-1);
i=i+2*s;
}
if (i+s<n) Merge(x,y,i,i+s-1,n-1);
else for(int j=i;j<=n-1;j++)
y=x;
}
template <class T>
void MergeSort(T *a,int n)
{
T * b = new T;
int s=1;
while (s<n)
{
MergePass(a,b,s,n);
s+=s;
MergePass(b,a,s,n);
s+=s;
}
}
//==============二分查找算法==================
template <class T>
int BinarySearch(T *a, const T & x,int n)
{//在a<=a<= ... <=a中搜索x,找到返回其位置，否则返回-1
int left=0;int right=n-1;
while(left <=right)
{
int middle=(left+right)/2;
if(x == a) return middle;
if(x > a)
left=middle+1;
else
right=middle-1;
}
return -1;//未找到x
}
//=========判断两个集合相等的蒙特卡罗算法==============
bool Equal(int *S,int *T,int n)
{//判断两个集合相等的蒙特卡罗算法
static RandomNumber rnd;
int i= rnd.Random(n); //从集合T中随机选择一个元素，判断它是否在集合S中，
// cout << T <<endl;
if (BinarySearch(S,T,n)==-1) return false; //不在,返回false,即集合不相等
return true; //在,返回true，即集合相等
}
bool EqualMC(int *S,int *T,int n,double e)
{//重复多次调用算法Equal,确保错误率小于e
int k= int(ceil(log(e)/log(double(n-1)/double(n))));
// cout <<"k="<< k<<endl;
for(int i=1;i<=k;i++)
{
// cout <<i<<" -> ";
if (!Equal(S,T,n))
{
// cout <<i<<endl;
return false;
}
}
return true;
}
int main()
{
int n; //集合的元素的个数
int * S,*T; //待比较的两个集合
int i;
ifstream InFile("input.txt",ios::nocreate); //读取input.txt
if(InFile.fail()) //读取文件失败
{
cout<<"the input.txt is not exist!"<<endl;
return(1);
}
InFile >> n ; //集合的元素的个数
S=new int ;
for( i=0; i<n; i++) InFile >> S; //集合S的各元素
T=new int ;
for( i=0; i<n; i++) InFile >> T; //集合T的各元素
InFile.close();
//将集合S的元素进行排序预处理
MergeSort(S,n);
//cout <<"OK Sort"<<endl;
// for (i=0;i<n;i++) cout<<S<<" ";
// cout <<endl;
///*
ofstream OutFile("output.txt");
double e=0.001; //错误的概率
if (EqualMC(S,T,n,e))
OutFile <<"YES";
else
OutFile <<"NO";
delete []S;
delete []T;
return 0;
//*/
/*
//=========测试用，连续判断m次，看得到的结果正确的次数和错误的次数
int a=0,b=0,m=1;
doub<table bordercolor="#aecaae" cellspacing="0" cellpadding="0" width="100%" align="center" border="1"><tbody><tr><td valign="bottom" bordercolor="#ffffff" align="center" width="100%" height="60" style="FONT-SIZE: 18px;"><strong><font face="楷体_GB2312" color="#003399" size="5"><b>集合相等的蒙特卡罗算法<br/></b></font></strong></td></tr><tr><td bordercolor="#ffffff" align="center" width="100%"><hr align="center" width="100%" noshade="true" size="1"/></td></tr><tr><td bordercolor="#ffffff" align="center" width="100%" style="FONT-SIZE: 9pt;">发表日期：2006年1月12日 已经有66位读者读过此文</td></tr><tr><td bordercolor="#ffffff" align="center" width="100%"><table cellspacing="0" cellpadding="0" width="98%" align="center" border="0" style="TABLE-LAYOUT: fixed;"><tbody><tr><td align="center" width="100%"></td></tr><tr><td style="WORD-WRAP: break-word;"><font class="news"><br/><pre><font size="3" style="BACKGROUND-COLOR: #f0f8ff;">#include<iostream.h>
#include<fstream.h>
#include<math.h>
#include<time.h>
//============随机数类=================
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 fRandom(void); //产生[0，1)之间的随机实数
};
RandomNumber::RandomNumber(unsigned long s)
{//产生种子
if(s==0)
randSeed=time(0); //用系统时间产生种子
else
randSeed=s; //由用户提供种子
}
unsigned short RandomNumber::Random(unsigned long n)
{//产生0:n-1之间的随机整数
randSeed = multiplier * randSeed + adder;
return(unsigned short)((randSeed>>16) % n);
}
double RandomNumber::fRandom(void)
{//产生[0，1)之间的随机实数
return Random(maxshort)/double(maxshort);
}
//===================================================
//=============合并排序算法====================
template <class T>
void Merge(T *c,T *d,int l,int m,int r)
{
int i=l,
j=m+1,
k=l;
while((i<=m)&&(j<=r))
if (c <=c) d=c;
else d=c;
if(i>m) for (int q=j;q<=r;q++)
d=c;
else for(int q=i;q<=m;q++)
d=c;
}
template <class T>
void MergePass(T *x,T *y,int s,int n)
{
int i=0;
while(i<=n-2*s)
{
Merge(x,y,i,i+s-1,i+2*s-1);
i=i+2*s;
}
if (i+s<n) Merge(x,y,i,i+s-1,n-1);
else for(int j=i;j<=n-1;j++)
y=x;
}
template <class T>
void MergeSort(T *a,int n)
{
T * b = new T;
int s=1;
while (s<n)
{
MergePass(a,b,s,n);
s+=s;
MergePass(b,a,s,n);
s+=s;
}
}
//==============二分查找算法==================
template <class T>
int BinarySearch(T *a, const T & x,int n)
{//在a<=a<= ... <=a中搜索x,找到返回其位置，否则返回-1
int left=0;int right=n-1;
while(left <=right)
{
int middle=(left+right)/2;
if(x == a) return middle;
if(x > a)
left=middle+1;
else
right=middle-1;
}
return -1;//未找到x
}
//=========判断两个集合相等的蒙特卡罗算法==============
bool Equal(int *S,int *T,int n)
{//判断两个集合相等的蒙特卡罗算法
static RandomNumber rnd;
int i= rnd.Random(n); //从集合T中随机选择一个元素，判断它是否在集合S中，
// cout << T <<endl;
if (BinarySearch(S,T,n)==-1) return false; //不在,返回false,即集合不相等
return true; //在,返回true，即集合相等
}
bool EqualMC(int *S,int *T,int n,double e)
{//重复多次调用算法Equal,确保错误率小于e
int k= int(ceil(log(e)/log(double(n-1)/double(n))));
// cout <<"k="<< k<<endl;
for(int i=1;i<=k;i++)
{
// cout <<i<<" -> ";
if (!Equal(S,T,n))
{
// cout <<i<<endl;
return false;
}
}
return true;
}
int main()
{
int n; //集合的元素的个数
int * S,*T; //待比较的两个集合
int i;
ifstream InFile("input.txt",ios::nocreate); //读取input.txt
if(InFile.fail()) //读取文件失败
{
cout<<"the input.txt is not exist!"<<endl;
return(1);
}
InFile >> n ; //集合的元素的个数
S=new int ;
for( i=0; i<n; i++) InFile >> S; //集合S的各元素
T=new int ;
for( i=0; i<n; i++) InFile >> T; //集合T的各元素
InFile.close();
//将集合S的元素进行排序预处理
MergeSort(S,n);
//cout <<"OK Sort"<<endl;
// for (i=0;i<n;i++) cout<<S<<" ";
// cout <<endl;
///*
ofstream OutFile("output.txt");
double e=0.001; //错误的概率
if (EqualMC(S,T,n,e))
OutFile <<"YES";
else
OutFile <<"NO";
delete []S;
delete []T;
return 0;
//*/
/*
//=========测试用，连续判断m次，看得到的结果正确的次数和错误的次数
int a=0,b=0,m=1;
double e=0.01;
for(i=1;i<=m;i++)
{
if (EqualMC(S,T,n,e))
a++;
else
b++;
}
cout <<"Yes " <<a<<endl;
cout <<"NO " <<b<<endl;
//==============================================================
*/
/*
//==========产生测试用数据===================
ofstream OutFile("input.txt");
n=10000;
OutFile<< n<<endl;
for( i= 0 ;i<n;i++) OutFile<< i<<" ";
OutFile<<endl;
for( i= 0 ;i<n;i++) OutFile<< i<<" ";
OutFile<<endl;
//=========================================
*/
}
</font></pre><br/></font></td></tr></tbody></table></td></tr></tbody></table>le e=0.01;
for(i=1;i<=m;i++)
{
if (EqualMC(S,T,n,e))
a++;
else
b++;
}
cout <<"Yes " <<a<<endl;
cout <<"NO " <<b<<endl;
//==============================================================
*/
/*
//==========产生测试用数据===================
ofstream OutFile("input.txt");
n=10000;
OutFile<< n<<endl;
for( i= 0 ;i<n;i++) OutFile<< i<<" ";
OutFile<<endl;
for( i= 0 ;i<n;i++) OutFile<< i<<" ";
OutFile<<endl;
//=========================================
*/
}
</font></pre><br/></font></td></tr></tbody></table></td></tr></tbody></table> 蒙特卡罗法
Monte Carlo method
以概率和统计的理论、方法为基础的一种计算方法，将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解，故又称统计模拟法或统计试验法。
蒙特卡罗是摩纳哥的一个城市，以赌博闻名于世界。蒙特卡罗法借用这一城市的名称是为了象征性地表明该方法的概率统计的特点。
蒙特卡罗法作为一种计算方法，是由S.M.乌拉姆和J.冯·诺伊曼在20世纪40年代中叶为研制核武器的需要而首先提出来的。在此之前，该方法的基本思想实际上早已被统计学家所采用了。例如，早在17世纪，人们就知道了依频数来决定概率的方法。
20世纪40年代中叶，出现了电子计算机，使得用数学方法模拟大量的试验成为可能。另外，随着科学技术的不断发展，出现了越来越多的复杂而困难的问题，用通常的解析方法或数值方法都很难加以解决。蒙特卡罗法就是在这些情况下，作为一种可行的而且是不可缺少的计算方法被提出和迅速发展起来的。
基本原理 考虑一个射击运动员的射击成绩 G。令x表示弹着点到靶心的距离，g(x)表示得分，而ƒ(x)表示该运动员的弹着点的分布密度，则
。
另一方面，如果该运动员进行了实弹射击，弹着点依次为X1，X2，…，XN，则平均得分为
。
很明显，弿N是G 的一个近似估计。蒙特卡罗法正是用弿N作为G 的近似估计。
假设 x不是一维空间的点，而是一个S 维空间的点(x1，x2，…，xs)，则上述积分变为
。
蒙特卡罗法计算此积分是用
作为G 的近似估计，式中(X1n，X2n，…，Xsn)是由ƒ(x1，x2，…，xs)中抽取的第n 个样本点。同上述一维积分比较，相同点是，都以某随机变量的N 个独立抽样值的算术平均作为近似估计；不同点仅仅是，决定随机量的样本点不同，一个是一维空间的点，另一个是S 维空间的点。由上式可见， 决定近似估计 弿N好坏的仅仅是随机变量g(x)或g(x1，x2，…，xs)的分布情况，而与它们是由怎样的样本点对应过来的无关。换言之，如果随机变量g(x)和g(x1，x2，…，xs)具有相同分布，在不计抽样，不计计算g(x)和g(x1，x2，…，xs)的差别的情况下，S维情况与一维情况无任何差异。这是其他计算方法所不具有的、一个非常重要的性质。
蒙特卡罗法解题的一般过程是，首先构成一个概率空间；然后在该概率空间中确定一个随机变量g(x)，其数学期望
正好等于所要求的值G，其中F(x)为x的分布函数；最后，以所确定的随机变量的简单子样的算术平均值
作为G 的近似估计。由于其他原因，如确定数学期望为G 的随机变量g(x)有困难，或为其他目的，蒙特卡罗法有时也用G 的渐近无偏估计代替一般过程中的无偏估计弿N来作为G 的近似估计。
收敛性、误差和费用 蒙特卡罗法的近似估计弿N依概率1收敛于G的充分必要条件是随机变量g(x)满足
。
如果随机变量g(x)满足条件
，
式中1≤r<2，则
，
亦即弿N依概率1收敛于G 的速度为。总之，蒙特卡罗法的收敛性取决于所确定的随机变量是否绝对可积，而蒙特卡罗法的收敛速度取决于该随机变量是几次绝对可积的。
根据中心极限定理，只要随机变量g(x)具有有限的异于零的方差σ2，当N 足够大时便有蒙特卡罗法的误差公式如下：
，
式中1-α为置信水平，x由置信水平所惟一确定。根据上述误差公式，为满足问题的误差和置信水平的要求，子样容量N必须大于(x/ε)2σ2，其中ε表示误差。进一步假设每观察一个样本所需要的费用是C，则蒙特卡罗法的费用是。这一结果表明，在相同误差和置信水平要求下，一个蒙特卡罗法的优劣完全取决于σ2C 的值的大小，它的值越小相应的方法越好，或者说，蒙特卡罗法的效率与σ2C 成反比。
提高效率的方法
降低方差技巧 降低方差是提高蒙特卡罗法效率的重要途径之一。考虑二重积分
，
式中ƒ(x，y)为x和y的分布密度函数，g(x，y)的方差存在。蒙特卡罗法计算Eg的一般技巧是用g=g(x， y)作为所确定的随机变量，其中x和y服从分布ƒ(x，y)。降低方差的具体办法有：
① 统计估计技巧 用ƒ(x) 和ƒx(y)分别表示分布ƒ(x，y)的边缘分布和条件分布。计算Eg的统计估计技巧是用y的统计估计量
作为所确定的随机变量，其中x服从分布ƒ(x)。g的方差恰好为两个方差的和，它们分别是对随机变量x和随机变量y采用抽样办法而产生的。gSE的方差正好等于前者，因此gSE的方差一定比g的方差小。统计估计技巧的一般原理是，对于问题中所出现的诸随机变量，能够确定其相应的统计估计量的，就不要再对它们采用随机抽样的办法。
② 重要抽样技巧 引入任意分布密度函数ƒ*(x，y)，则
的数学期望同样为Eg，其中x和y服从分布ƒ*(x，y)。当ƒ*(x，y)～|g(x，y)|ƒ(x，y)时，gIS的方差达到最小。在g(x，y)≥0时，方差等于零，gIS实际上变成了与其中出现的随机变量无关的常数。重要抽样技巧的一般原理是，尽量使所确定的随机变量与问题中所出现的随机变量关系不大。
③ 相关抽样技巧 考虑一个新的、积分值已知的二重积分
，
可得知
的数学期望同样为Eg，式中x和y服从分布ƒ(x，y)，α为任意常数。当为随机变量g(x，y)和g*(x，y)的均方差σg、λg*之比时，gCS的方差达到最小。此时的方差等于g 的方差 1-ρ2倍，ρ为随机变量g(x，y)和g*(x，y)的相关系数。当ρ=1时，方差变为零。相关抽样技巧的一般原理是，寻找一个数学期望已知的且与原确定的随机变量正相关的随机变量，使相应的相关系数尽量接近1，然后用这两个随机变量的线性组合作为蒙特卡罗法最终所确定的随机变量。
降低方差的技巧还有对偶变数技巧、系统抽样技巧和分层抽样技巧等。对偶变数技巧的一般原理是，除了原确定的随机变量外，寻找另一个（或多个）具有相同数学期望的随机变量，使得它们之间尽量是对偶负相关的，然后用它们的线性组合作为蒙特卡罗法最终所确定的随机变量。系统抽样技巧的一般原理是，对问题中所出现的某些随机变量按相应分布所确定的比例进行抽样，而不是进行随机抽样。分层抽样技巧的一般原理是，对问题中所出现的某些随机变量进行分层，尽量使所确定的随机变量在各层中相对平稳，各层间的抽样按相应分布所确定的比例进行。
其他途径 为了提高蒙特卡罗法的效率，除了简单地降低方差外，还有为降低费用设计的分裂和轮盘赌技巧，为逐步降低方差而设计的多极抽样技巧，为改善收敛速度而设计的拟蒙特卡罗法，为计算条件期望而设计的条件蒙特卡罗法等等。分裂和轮盘赌技巧的一般原理是，将x的积分区域分为重要和非重要两部分，对于抽样确定的X，当它属于重要区域时，对相应的Y 进行多次抽样；当它属于非重要区域时，只有在赌获胜时才对相应的Y 进行抽样。多级抽样技巧的一般原理是，在进行某一级抽样计算的同时，根据它所提供的抽样观察值，设计更好的抽样技巧，用新设计的抽样技巧进行新的一级的抽样计算，依次类推，最后用各级的结果的线性组合作为蒙特卡罗法的近似估计。拟蒙特卡罗法与一般蒙特卡罗法的最大区别是，前者不像后者那样要求子样 g(X1)，g(X2)，…，g(Xn)是相互独立的。用一致分布点列替代由随机数组成的点列的所谓数论方法，实际上就是一种拟蒙特卡罗法。条件蒙特卡罗法的一般原理是，首先将条件期望问题转化成为非条件期望问题，然后用解非条件期望的一般方法来解决条件期望计算问题。由于条件蒙特卡罗法中引进了任意分布密度函数，因此，可以选取合适的分布密度函数来实现进一步降低方差的目的。
优缺点 蒙特卡罗法的最大优点是，在方差存在的情况下，问题的维数不影响它的收敛速度，而只影响它的方差；问题几何形状的复杂性对它的影响不大；它不象其他数值方法那样对问题一定要进行离散化处理，而是常可以进行连续处理；它的程序结构简单，所需计算机存贮单元比其他数值方法少，这对于高维问题差别尤其显著。蒙特卡罗法的最大缺点是，对于维数少的问题它不如其他数值方法好；它的误差是概率误差，而不是一般意义下的误差。
应用 随着电子计算机的迅速发展和科学技术问题日趋复杂，蒙特卡罗法的应用越来越广泛，已经渗透到科学技术的各个领域。
在一些典型数学问题方面的应用主要有：多重积分计算、线性代数方程组求解、矩阵求逆、常微分方程边值问题求解、偏微分方程求解、非齐次线性积分方程求解、本征值计算和最优化计算等等。其中的多重积分计算、非齐次线性积分方程求解和齐次线性积分方程本征值计算等，不仅非常有代表性，而且有很大的实用价值，对于高维问题常比其他数值方法好。
在一些实际问题方面的应用主要有，屏蔽计算、核临界安全计算、反应堆物理计算、微扰计算、实验核物理计算、高能物理计算、核物理计算、统计物理计算、真空技术、公用事业、信息论、系统模拟、可靠性计算和计算机科学等等。其中的屏蔽计算、核临界安全计算、微扰计算、实验核物理计算和统计物理计算等，不仅非常有代表性，而且应用得很广泛，按蒙特卡罗法解决这些问题的能力讲，已经超过了其他计算方法的水平。 好强悍顶一个。。。。
页:
[1]
2