蒙特卡洛理解与python实现2

蒙特卡洛定义

模拟是指把某一现实的或抽象的系统的某种特征或部分状态， 用另一系统（称为模拟模型）来代替或模拟。 为了解决某问题， 把它变成一个概率模型的求解问题， 然后产生符合模型的大量随机数， 对产生的随机数进行分析从而求解问题， 这种方法叫做随机模拟方法， 又称为蒙特卡洛(Monte Carlo)方法。

蒙特卡洛方法的名字来源于摩纳哥的一个城市蒙特卡洛，该城市以赌博业闻名，而蒙特卡洛方法正是以概率为基础的方法，第一次提出是在二战时期曼哈顿计划制造第一代原子弹的项目中。蒙特卡洛方法是一系列符合该定义算法的统称，而不是一个算法。

蒙特卡洛主要有三个用处：

• 近似求解积分（蒙特卡洛积分法）
• 密度估计：通过采样逼近目标分布
• 优化函数：找到能使目标函数最大化或最小化的样本

初识蒙特卡洛- 随机投点法计算圆周率

投点法是最朴素的蒙特卡洛思想，与几何概率类似。

image.png

算法实践：

def estimate_pi():
    N = 10000
    count = 0
    for i in range(N):
        x = np.random.uniform(0,1)
        y = np.random.uniform(0,1)
        if x**2+y**2 <= 1:
            count += 1
    return (count/N)*4

标准蒙特卡洛-平均值法

英文称呼有以下几种：standard monte carlo，directed monte carlo， standard sample average 。
原理是利用强大数定理，n足够大时，样本的均值等于随机变量的期望。

求积分

假设p(x)是目标分布，我们要计算, 则使用标准蒙特卡洛法估算期望

得

其中,按照强大数定理，依概率为1收敛到。
通常

image.png

• 算法实践：利用平均采样法求

1. 生成随机变量U
2. 采样N次 代入

def average_sample():
    N = 1000
    x_hat = []
    for i in range(N):
        u = np.random.uniform(0,1)
        x_hat.append(np.exp(u))
    return sum(x_hat)/N

重要性采样

使用标准蒙特卡洛法有时候生成的结果方差比较大，重要性抽样能够有效的减小方差。
标准的蒙特卡洛法是
加上一个重要性分布,重写以上公式得

我们可以得到新的估计

其中
根据强大数定理SLLN，依概率为1收敛到

self-normlized重要性采样

有时不仅很难直接抽样，而且本身未知， 只能确定到差一个常数倍的，计算c一般很困难。定义,则得到新的蒙特卡洛积分公式

如何选择一个好的重要性分布q- Tilted Densities方法

Tilted Densities方法利用moment generating function ，构造如下：

• 算法实践 利用重要性采样法求

1. 对被积函数p(x)做泰勒展开得
2. 取重要性分布
3. 利用逆变换法产生
   q(x)的分布函数Q(x)求逆得
   生成随机变量,令
4. 则重要抽样法的积分公式为
   

def importance_sample():
    N = 1000
    x_hat = []
    for i in range(N):
        x = np.random.uniform(0,1) #以0-1均匀分布生成随机数x
        x = np.sqrt(1+3*x) - 1 #将x代入，生成PDF为q(x)的分布
        y = np.exp(x) / ((2/3)*(1+x)) #以p(x)/q(x)采样
        x_hat.append(y) 
    return sum(x_hat)/N   #对所有仿真结果求均值 

• 算法实践2：假设, 估计概率的值

引入指示变量I,则有

1. 使用标准蒙特卡洛方法
   取N个样本 则
2. 使用重要性抽样方法
   其中, ，我们选择,
   则
   其中,当时候q(x)是一个合理的选择

def average_sample():
    count = 0
    for i in range(1000):
        x_k = np.random.normal(0,1)
        if x_k > 8:
            count += 1
    return count/1000

def importance_sample():
    N = 50000
    alpha = 8
    result = []
    for i in range(N):
        x = np.random.normal(alpha,1) #q分布 N(alpha,1)
        if x > alpha:
            result.append(np.exp(0.5*alpha**2 - alpha*x))
    return sum(result)/N

参考文献

1. http://www.columbia.edu/~mh2078/MonteCarlo.html关于重要性采样这里有更数学化的证明过程以及其他选择q(x)的方法。

2. 北京大学统计计算讲义-李东风 本文中部分例题来源于此，里面提供了关于重要性采样如何减小方差的理论验证过程。

3. https://people.duke.edu/~ccc14/sta-663/MonteCarlo.html
   这里提供了蒙特卡洛方法的纯numpy运算实现，不用for循环，代码更简洁，但没那么好理解。