第十课.蒙特卡洛方法与随机过程

蒙特卡洛方法

蒙特卡洛的概念

蒙特卡洛是一个赌场的名称，用它作为名字是因为蒙特卡洛方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡洛方法都是为了求解一些不太好求解的求和或者积分问题。

例如下图是一个经典的用蒙特卡洛求圆周率的问题，用计算机在一个正方形之中随机的生成点，计数有多少点落在$1/4$圆之中，这些点的数目除以总的点数目即圆的面积，根据圆面积公式即可求得圆周率：

蒙特卡洛方法的另一个应用是求积分，某些函数的积分不好求，我们可以按照下面的方法将这个函数进行分解，之后转化为求期望与求均值的问题：
$${\int }_a^{b}h(x)dx={\int }_a^{b}f(x)p(x)dx=E_{p(x)}[f(x)]$$
从分布$p(x)$中采样大量的样本$x_1,x_2,x_3,...,x_N$，这些样本符合分布$p(x)$，于是代入样本得到：
$${\int }_a^{b}h(x)dx=E_{p(x)}[f(x)]=\frac{1}{N}\sum _i^{N}f(x_i)$$

蒙特卡洛采样方法

对于概率分布$p(x)$的采样又分为3种：

• 直接采样；
• 接受-拒绝采样；
• 重要性采样。

直接采样

直接采样的方法是根据概率分布进行采样。对一个已知概率密度函数与累积概率密度函数的概率分布，我们可以直接从累积分布函数（cdf）进行采样。

如下图所示是四个高斯分布的累积概率分布函数，可以看出函数的值域是$(0,1)$，我们可以从$U(0,1)$均匀分布中进行$N$次采样，再根据某个高斯分布的累积分布函数的反函数计算出对应的$x$，这样就获得了符合某个高斯分布的$N$个粒子。

使用累积分布函数进行采样看似简单，但是由于很多分布我们并不能写出概率密度函数与累积分布函数，所以这种方法的适用范围较窄。

接受-拒绝采样

对于累积分布函数未知的分布，我们可以采用接受-拒绝采样。如下图所示，$p(z)$是我们希望采样的分布，$q(z)$是我们提议的分布(proposal distribution)，令$kq(z)>p(z)$，我们首先在$kq(z)$中按照直接采样的方法采样粒子，接下来判断这个粒子落在途中什么区域，对于落在灰色区域的粒子予以拒绝，落在红线下的粒子接受，最终得到符合$p(z)$的$N$个粒子。

重要性采样

接受拒绝采样完美的解决了累积分布函数不可求时的采样问题。但是接受拒绝采样非常依赖于提议分布(proposal distribution)的选择，如果提议分布选择的不好，可能采样时间很长却获得很少满足分布的粒子。而重要性采样就解决了这一问题。

直接采样与接受拒绝采样都是假设每个粒子的权重相等，而重要性采样则是给予每个粒子不同的权重，使用加权平均的方法来计算期望：
$$E_{p(x)}[f(x)]={\int }_a^{b}f(x)\frac{p(x)}{q(x)}q(x)dx=E_{q(x)}[f(x)\frac{p(x)}{q(x)}]$$
我们从提议分布$q(x)$中采样粒子$x_1,...,x_N$，每个粒子的权重为$p(x_i)/q(x_i)$，从而得到：
$$E_{p(x)}[f(x)]=\frac{1}{N}\sum _i^{N}f(x_i)\frac{p(x_i)}{q(x_i)}$$

小结

蒙特卡洛方法是一种近似推断的方法，通过采样大量粒子的方法来求解期望、均值、面积、积分等问题，蒙特卡洛对某一种分布的采样方法有直接采样、接受拒绝采样与重要性采样三种，直接采样最简单，但是需要已知累积分布的形式。接受拒绝采样与重要性采样适用于原分布未知的情况，这两种方法都是给出一个提议分布，不同的是接受拒绝采样对不满足原分布的粒子予以拒绝，而重要性采样则是给予每个粒子不同的权重。

随机过程

随机变量与随机过程

随机过程是一串随机变量的序列，在这个序列中，每一个数据都是一个随机变量，因此在随机过程的概率模型处理中，重点关注时间和数据两方面内容，所以随机过程就是一串（有限或者无限的）随机变量序列。常见的随机过程建模场景无数，比如：

• 足球赛中每场比赛进球数构成的序列；
• 某个路口每分钟通过的车辆数量构成的序列。

在讨论随机变量时，关注的是其分布；讨论随机过程时，将随机变量连成序列，关注的是序列之间的相关关系，序列整体的统计特征。

随机过程实际场景举例

下面举例两个实际场景，用随机过程对其建模，使用蒙特卡洛方法对过程进行模拟。

赌博中的随机过程

赌徒和商人对赌抛硬币，如果为正面朝上，则本轮赌徒获胜，商人付给赌徒1美元，如果反面朝上，则本轮赌徒失败，付给商人1美元。赌徒有10美元的资本，手上的钱一旦为0，就退出赌局。

分析这个赌博的过程，赌徒的赌博结果本质上取决于每次投掷硬币的结果，每一轮赌博就是一个伯努利实验，获胜的概率为0.5，赌博的过程就是一串伯努利实验构成的随机过程，每轮赌局，获胜则资本加1，失败则资本减1。

对于一个赌徒，一旦进入赌局，每轮赌局的结果构成的序列是唯一的，我们想要看到赌博过程背后的特征，于是要用蒙特卡洛方法进行模拟，使用大量样本，观察整体特征。比如设置赌徒数量为10，每个赌徒最多赌博100轮：如果提前输光资本，就必须退出，如果到第100轮还有资本，也要停止。

模拟过程如下：

import pandas as pd
import random

sample_list=[]
round_num=100
person_num=10

for person in range(1,person_num+1):
    money=10
    for the_round in range(1,round_num+1):
    	# 等概率生成硬币正反面结果
        result=random.randint(0,1)
        if result==1:
            money=money+1
        elif result==0:
            money=money-1
        if money==0:
            break
    sample_list.append([person,the_round,money])
    
sample_df=pd.DataFrame(sample_list,columns=['person','round','money'])
sample_df.set_index('person',inplace=True)
print(sample_df)

结果为：

        round  money
person              
1         100     10
2         100      6
3         100     26
4          24      0
5         100     12
6         100      8
7         100     10
8          58      0
9         100     18
10         60      0

从中看出，10个人有6个提前输光，剩下的4人只有3人是盈利的。输赢的概率本来是一半，但赌徒破产走人的结局却时常发生，其原因在于商人的资本是无限的，而赌徒的资本有限。

股价变化的过程

再以股票为例，在金融行业，有一个公式：
$$S_{t+1}=S_t+\hat{\mu }{S}_t\Delta t+\sigma S_tϵ\sqrt{\Delta t}$$
公式利用目前的股价$S_t$去预测$\Delta t$时间之后的股价$S_{t+1}$，其中，$\hat{\mu }$表示股票收益率的期望，一般为0.15，$\sigma$为股票的波动率，一般为0.2，$ϵ$为服从标准正态分布的随机变量，$ϵ$导致股价$S_i$也成为随机变量，而股价构成的序列就是一个随机过程。

同样用蒙特卡洛方法，估计在$S_0=10$的情况下，244个$\Delta t$后的股价概率分布。模拟如下：

%matplotlib inline
import scipy
import matplotlib.pyplot as plt
from math import sqrt

s0=10.0
T=1.0
n=224*T
mu=0.15
sigma=0.2

n_simulation=10000

dt=T/n
s_array=[]

for i in range(n_simulation):
    s=s0
    for j in range(int(n)):
        e=scipy.random.normal()
        s=s+mu*s*dt+sigma*s*e*sqrt(dt)
    s_array.append(s)
    
plt.hist(s_array,bins=30,normed=True,edgecolor='k')
plt.show()

上面模拟了10000个样本经过上述随机过程后的股价分布情况。但这只能体现每个随机过程的最终股价，有时人们想要了解股价的变化过程，所以要对每个随机过程都进行记录：

%matplotlib inline
import scipy
import matplotlib.pyplot as plt
from math import sqrt
import numpy as np

s0=10.0
T=1.0
n=224*T
mu=0.15
sigma=0.2

n_simulation=100

dt=T/n
random_series=np.zeros(int(n),dtype=float)
x=range(0,int(n))

for i in range(n_simulation):
    random_series[0]=s0
    for j in range(1,int(n)):
        e=scipy.random.normal()
        random_series[j]=random_series[j-1]+mu*random_series[j-1]*dt+sigma*random_series[j-1]*e*sqrt(dt)
    plt.plot(x,random_series)
    
plt.show()

以上模拟了100个随机过程，每个随机过程的变化展现，有利于人们更合理地分析股价。

两个重要随机过程

随机过程种类繁多，应用广泛，人们常用的随机过程有两类：到达过程和马尔科夫过程；

到达过程：关注的是某种"到达"事件是否发生，比如：在一个服务窗口前，顾客的到达时刻，十字路口车辆依次通过的时刻。重点是：相邻两次到达时刻之间的时间（相邻间隔时间）为相互独立的随机变量。

当到达时间为离散情况时，随机过程为伯努利过程，相邻间隔时间服从几何分布。当到达时间为连续的情况，随机过程为泊松过程，相邻间隔时间服从指数分布。

这类过程的特点在于相邻间隔时间相互独立，随机过程无记忆性，比如赌博过程中，本轮赌局的输赢不对下一轮赌局的输赢带来影响。

---

几何分布：几何分布（Geometric distribution）是离散型概率分布。其中一种定义为：在$n$次伯努利试验中，试验$k$次才得到第一次成功的机率。详细地说，是：前$k-1$次皆失败，第$k$次成功的概率。

在伯努利试验中，记每次试验中事件A发生的概率为$p$，试验进行到事件A出现时停止，此时所进行的试验次数为$X$，其分布列（离散的概率密度）为：
$$p(X=k)=(1-p{)}^{k-1}p,k=1,2,...$$
记$X\sim GE(p)$。

指数分布：指数分布（也称为负指数分布）是描述泊松过程中的事件之间的时间的概率分布，即事件以恒定平均速率连续且独立地发生的过程。其概率密度函数为：
$$p(X=x)=\lambda e^{-\lambda x},x>0,p(X=x)=0,x\le 0$$
记作$X\sim E(\lambda )$，其中，$\lambda >0$是分布的率参数，即每单位时间内发生某事件的次数。$X$表示相邻间隔时间。

指数分布的期望为$1/\lambda$。

泊松分布：Poisson分布，是一种统计与概率学里常见到的离散概率分布， 泊松分布适合于描述单位时间内随机事件发生的次数。其概率密度函数为：
$$p(X=k)=\frac{{\lambda }^k}{k!}{e}^{-\lambda }$$
其中，$X$表示单位时间内随机事件的发生次数，泊松分布的期望为$\lambda$。从期望上看，泊松分布与指数分布正好呈倒数关系。

泊松过程是一个计数的过程，在$(0,t)$时与$(0,t+s)$时事件发生的次数独立且同服从参数为$\lambda$的泊松分布。这样按照时间进行的事件发生的次数就是一个泊松过程。

---

马尔科夫过程：与到达过程不同，未来的数据依赖当前的数据，比如股价变化过程，未来的股价与现在的股价有关，但注意，未来的数据仅仅与当前数据有关，与过去的历史数据无关。