蒙特卡罗算法

蒙特卡罗算法

定义

• • 
    
• 蒙特卡罗是一类随机方法的统称。这类方法的特点是，可以在随机采样上计算得到近似结果，随着采样的增多，得到的结果是正确结果的概率逐渐加大，但在（放弃随机采样，而采用类似全采样这样的确定性方法）获得真正的结果之前，无法知道目前得到的结果是不是真正的结果。
  • 
    

特点

• • 
    
• 基于随机数
• 统计模拟法
• 拟蒙特卡洛算法在近几年也获得迅速发展。这种方法是用确定性的超均匀分布代替蒙特卡洛算法中的
  随机数序列，对于某些特定问题计算速度比普通的蒙特卡洛算法高几百倍。
  • 
    

对比

• • 
    
• 拉斯维加斯算法
  尽量找最好的，但不保证能找到
• 蒙特卡罗算法
  尽量找好的，但不保证是最好的
  • 
    

实现的三个步骤

1. 1. 
      
2. 构造描述概率过程
3. 实现从已知概率分布抽样
   随机数是实现模拟的工具
4. 建立各种估计量
   得到无偏估计
   1. 
      

举例

• • 
    
• 模拟股票
• 
  

​​​https://zhuanlan.zhihu.com/p/28133690​​​

• 
  

clear
dt=1/365.0;                   % 一天的年单位时间
S0=20;                       % 股票在初始时刻的价格，程序中假设
r=0.031;                      % 期望收益率
sigma=0.6;                    % 波动率=0.6
expTerm=r*dt;                 % 漂移项dt
stddev=sigma*sqrt(dt);       % 波动项o:dz(t) 
nDays1=90;                    % 要模拟的总天数
for nDays=1:nDays1            % nDays表示时刻t
nTrials=10000;                % 模拟次数
for j=1:nTrials 
n = randn(1,nDays);           %生成nDays个标准正态分布随机数
S=S0; 
for i=1:nDays 
dS = S*(expTerm+stddev*n(i));   % 模拟计算股票价格的增量
S=S+dS;                      %计算股票价格
end
S1(nDays,j)=S;               % 将每天的股票模拟价格数据记录在S1中
end
end
S2=mean(S1');                 % 计算每天模拟的股票价格的均值，作为价格的估值 
plot(S2','-o')                    % 90天期间股票价格估值的曲线图
figure(2)
hist(S1(90,:),0:0.5:65)      %第90天的股票价格模拟的直方图

• • 
    
• 计算圆周率
  • 
    

#include <bits/stdc++.h>

#define MAX_ITERS 1000000

using namespace std;

double Rand(double L, double R)
{
    return L + (R - L) * rand() * 1.0 / RAND_MAX;
}

double GetPi()
{
    srand(time(NULL));
    int cnt = 0;
    for(int i = 0; i < MAX_ITERS; i++)
    {
        double x = Rand(-1, 1);
        double y = Rand(-1, 1);
        if(x * x + y * y <= 1)
            cnt++;
    }
    return cnt * 4.0 / MAX_ITERS;
}

int main()
{
    for(int i = 0; i < 10; i++)
        cout << GetPi() << endl;
    return 0;
}

• • 
    
• 计算积分
• 详细过程见：
  • 
    

#include <bits/stdc++.h>

#define MAX_ITERS 10000000

using namespace std;

struct Point
{
    double x, y;
};

double Rand(double L, double R)  
{  
    return L + (R - L) * rand() * 1.0 / RAND_MAX;  
} 

Point getPoint()
{
    Point t;
    t.x = Rand(1.0, 2.0);
    t.y = Rand(0.0, 1.0);
    return t;
}

double getResult()
{
    int m = 0;
    int n = MAX_ITERS;
    srand(time(NULL));
    for(int i = 0; i < n; i++)
    {
        Point t = getPoint();
        double res = t.x * t.y;
        if(res <= 1.0)
            m++;
    }
    return pow(2.0, 1.0 * n / m);
}

int main()
{
    for(int i = 0; i < 20; i++)
        cout << fixed << setprecision(6) << getResult() << endl;
    return 0;
}

代码改变世界