蒙特卡洛的个人整理

结论

结论先行：蒙特卡洛是一类随机算法，可以抽象为两步：①从分布中采样②后续操作。下面就详细看看这些例子。
注：本文主要按照自己的思路对知识进行串联，包含大量引用（见文末）

0. 常见的两种蒙特卡洛方式

0.1 投点法

从连续均匀分布中抽样，然后通过一个“操作”得到了落在面积内的点，然后依照古典概型计算概率。
这方式就是抽样+投点

import math
import random
m = 10000
n = 0
for i in range(m):
	# x、y为0-1之间的随机数
    x = random.random()
    y = random.random()
    # 若点(x,y) 属于图中1/4圆内 则有效个数+1
    if math.sqrt(x**2 + y**2) < 1:
        n += 1
# 计算pi
pi = 4 * n / m
print("pi = {}".format(pi))

# pi = 3.1508（结果具有随机性 不一定完全一样）

0.2 均值法

如果f(x)难以计算原函数，可以通过这样求均值的方式，来近似原来的积分。p(x)为自己随便假设的分布。这种类型就是：抽样+均值

1.从分布中采样

从分布中采样可以从两个角度来看：从常规分布中采样，比如：正态分布、均匀分布等或者从不常规的分布中采样，比如混合模型。

1.1 从常规分布中采样

1.2 从不常规的分布中采样（接受拒绝采样）

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, uniform

#目标采样分布的概率密度函数
def p(x):
    return 0.3*np.exp(-(x-0.3)**2)+0.7*np.exp(-(x-2.)**2/0.3)

#建议分布G：选择均值为1.4，方差为1.2的正态分布
G_rv = norm(loc=1.4, scale=1.2)

#常数值C=2.5
C=2.5

#均匀分布U(0,1)
uniform_rv = uniform(loc=0, scale=1)

#绘制目标分布的概率密度函数p(x)与建议分布G的概率密度函数g(x)的C倍
x = np.arange(-4., 6., 0.01)
plt.plot(x, p(x), color='r', lw=2, label='p(x)')
plt.plot(x, C*G_rv.pdf(x), color='b', lw=2, label='C*g(x)')  #g.pdf(x)表示正态分布的概率密度函数
plt.legend()
plt.show()

sample = []
#设10000个候选采样点
for i in range(10000):
    #step1：从建议分布G中进行采样，得到服从G分布的随机数
    Y = G_rv.rvs(1)[0]   #rvs():产生服从指定分布的随机数
    
    #step2：从均匀分布 U(0,1) 中进行采样，得到服从标准均匀分布的随机数
    U = uniform_rv.rvs(1)[0] 
    
    #step3：判断，如果 P(Y)≥U*C*g(Y)，则接受
    if p(Y)>=U*C*G_rv.pdf(Y):
        sample.append(Y)


#绘制目标分布的概率密度函数p(x)与建议分布G的概率密度函数g(x)的C倍
x = np.arange(-3., 5., 0.01)
plt.gca().axes.set_xlim(-3, 5)
plt.plot(x, p(x), color='r')
plt.hist(sample, color='b', bins=200, density=True, stacked=True, edgecolor='b')  #把所有直方归一化到1
plt.show()

参考文献

1.https://zhuanlan.zhihu.com/p/369099011
2.https://haihong.blog.csdn.net/article/details/119932032?spm=1001.2101.3001.6650.3&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-3-119932032-blog-79825783.235%5Ev32%5Epc_relevant_default_base3&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromBaidu%7ERate-3-119932032-blog-79825783.235%5Ev32%5Epc_relevant_default_base3&utm_relevant_index=4