# 强化学习-学习笔记6 | 蒙特卡洛算法

Monte Carlo Algorithms. 蒙特卡洛算法是一大类随机算法，又称为随机抽样或统计试验方法，通过随机样本估计真实值。

## 6. 蒙特卡洛算法

### 6.1 计算 $$\pi$$

#### a. 原理

$$\frac{4m}{n}\rightarrow \pi$$，as n → ∞

$$|\frac{4m}{n}-\pi|=O(\frac{1}{\sqrt{n}})$$

#### b. 代码

#coding=utf-8
#蒙特卡罗方法计算 pi
import random,math,time
start_time = time.perf_counter()
s = 1000*1000
hits = 0
for i in range(s):
x = random.random()
y = random.random()
z = math.sqrt(x**2+y**2)
if z<=1:
hits +=1

PI = 4*(hits/s)
print(PI)
end_time = time.perf_counter()
print("{:.2f}S".format(end_time-start_time))

# 输出
3.141212
0.89S


### 6.2 Buffon's Needle Problem

#### a. 原理

6.1 类似，我们随机扔 n 根针，这样相交个数的期望$$Pn = \frac{2ln}{\pi{d}}$$ 。我们可以观察到（如果是电脑模拟即为通过公式判断出）有 m 跟针实际与线相交，如果n足够大，则 $$m\approx \frac{2ln}{\pi{d}}$$

$$\pi$$ 公式即为： $$\pi\approx \frac{2ln}{md}$$

#### b. 代码

import numpy as np

def buffon(a,l,n):
xl = np.pi*np.random.random(n)
yl = 0.5*a*np.random.random(n)
m = 0
for x,y in zip(xl,yl):
if y < 0.5*l*np.sin(x):
m+=1
result = 2*l/a*n/m
print(f'pi的估计值是{result}')

buffon(2,1,1000000)

# 输出为：
pi的估计值是3.153977165205324


import matplotlib.pyplot as plt
import random
import math
import numpy as np

NUMBER_OF_NEEDLES = 5000

class DefineNeedle:
def __init__(self, x=None, y=None, theta=None, length=0.5):
if x is None:
x = random.uniform(0, 1)
if y is None:
y = random.uniform(0, 1)
if theta is None:
theta = random.uniform(0, math.pi)

self.needle_coordinates = np.array([x, y])
self.complex_representation = np.array(
[length/2 * math.cos(theta), length/2*math.sin(theta)])

def intersects_with_y(self, y):
return self.end_points[0][1] < y and self.end_points[1][1] > y

class BuffonSimulation:
def __init__(self):
self.floor = []
self.boards = 2
self.list_of_needle_objects = []
self.number_of_intersections = 0

fig = plt.figure(figsize=(10, 10))
self.buffon = plt.subplot()
self.results_text = fig.text(
0, 0, self.estimate_pi(), size=15)
self.buffon.set_xlim(-0.1, 1.1)
self.buffon.set_ylim(-0.1, 1.1)

def plot_floor_boards(self):
for j in range(self.boards):
self.floor.append(0+j)
self.buffon.hlines(
y=self.floor[j], xmin=0, xmax=1, color='black', linestyle='--', linewidth=2.0)

def toss_needles(self):
needle_object = DefineNeedle()
self.list_of_needle_objects.append(needle_object)
x_coordinates = [needle_object.end_points[0]
[0], needle_object.end_points[1][0]]
y_coordinates = [needle_object.end_points[0]
[1], needle_object.end_points[1][1]]

for board in range(self.boards):
if needle_object.intersects_with_y(self.floor[board]):
self.number_of_intersections += 1
self.buffon.plot(x_coordinates, y_coordinates,
color='green', linewidth=1)
return
self.buffon.plot(x_coordinates, y_coordinates,
color='red', linewidth=1)

def estimate_pi(self, needles_tossed=0):
if self.number_of_intersections == 0:
estimated_pi = 0
else:
estimated_pi = (needles_tossed) / \
(1 * self.number_of_intersections)
error = abs(((math.pi - estimated_pi)/math.pi)*100)
return (" Intersections:" + str(self.number_of_intersections) +
"\n Total Needles: " + str(needles_tossed) +
"\n Approximation of Pi: " + str(estimated_pi) +
"\n Error: " + str(error) + "%")

def plot_needles(self):
for needle in range(NUMBER_OF_NEEDLES):
self.toss_needles()
self.results_text.set_text(self.estimate_pi(needle+1))
if (needle+1) % 200 == 0:
plt.pause(1/200)
plt.title("Estimation of Pi using Probability")

def plot(self):
self.plot_floor_boards()
self.plot_needles()
plt.show()

simulation = BuffonSimulation()
simulation.plot()


### 6.3 估计阴影部分的面积

$\begin{cases} x^2+y^2>4\\ (x-1)^2+(y-1)^2\leq1\end{cases}$

• 易知，正方形面积 $$A_1=4$$；设阴影部分面积为 $$A_2$$
• 随机抽样的点落在阴影部分的概率为：$$P=\frac{A_2}{A_1}=\frac{A_2}{4}$$
• 从正方形区域抽样 n 个点，n尽可能大，则来自阴影部分点的期望为：$$nP=\frac{nA_2}{4}$$
• 如果实际上满足上述条件的点 有 m 个，则令 $$m\approx nP$$
• 得到：$$A_2\approx \frac{4m}{n}$$

### 6.4 求不规则积分

• 从区间 $$[a,b]$$ 上随机均匀抽样 $$x_1,x_2,...,x_n$$;

• 计算 $$Q_n = (b-a)\frac{1}{n}\sum_{i=1}^nf(x_i)$$，即均值乘以区间长度；

这里均值乘以区间长度是 实际值，而 I 是期望值

• $$Q_n$$ 近似 $$I$$

• 从区间 $$\Omega$$ 上随机均匀抽样 $$\vec{x_1},\vec{x_2},...,\vec{x_n}$$;

• 计算 $$\Omega$$ 的体积V（高于三维同样）：$$V=\int_\Omega{d\vec{x}}$$

hh值得注意的是，这一步仍要计算定积分，如果形状过于复杂，无法求得 V，那么无法继续进行，则无法使用蒙特卡洛算法。所以只能适用于比较规则的区域，比如圆形，长方体等。

• 计算 $$Q_n =V \frac{1}{n}\sum_{i=1}^nf(\vec{x_i})$$，即均值乘以区间长度；

这里均值乘以区间长度是 实际值，而 I 是期望值

• $$Q_n$$ 近似 $$I$$

• 定义一个二元函数 $$f(x,y)=\begin{cases} 1 \ \ if点在圆内\\ 0 \ \ if 点在圆外\end{cases}$$;
• 定义一个区间 $$\Omega=[-1,1]×[-1,1]$$
• $$I =\pi {r^2}=\pi$$
• 接下来用蒙特卡洛近似 I，得到关于 $$\pi$$的算式即可得到近似的$$\pi$$;
• 随机抽样 n 个点，记为$$(x_1,y_1),...,(x_n,y_n)$$
• 计算 区域面积 $$V = \int_\Omega{dxdy}=4$$
• 计算 $$Q_n =V \frac{1}{n}\sum_{i=1}^nf(x_i,y_i)$$
• 蒙特卡洛近似 Q 与 I 近似相等：$$\pi=Q_n=\int_\Omega{f(x,y)}{dxdy}$$

### 6.5 用蒙特卡洛近似期望

• 定义 X 是 d 维的随机变量，函数 p(x) 是一个PDF，概率密度函数；
• 函数 $$f(x)$$ 的期望：$$\mathbb{E}_{x\sim{p}}[f(X)]=\int_{R^d}f(X)\cdotp(x)dx$$
• 直接以上面的方式求期望可能并不容易，所以通常使用蒙特卡洛近似求期望：
1. 随机抽样：根据概率密度函数 $$p(x)$$ 进行随机抽样，记为$$X_1,X_2,...,X_n$$
2. 计算 $$Q_n =\frac{1}{n}\sum_{i=1}^nf(x_i)$$
3. 用 Q 近似 期望$$\mathbb{E}_{x\sim{p}}[f(X)]$$

### 6.6 总结 | 蒙特卡洛算法的思想

• 蒙特卡洛是摩洛哥的赌场；

• 蒙特卡洛算法得到的结果通常是错误的，但很接近真实值，对于对精度要求不高的机器学习已经足够。

随机梯度下降就是一种蒙特卡洛算法，用随机的梯度近似真实的梯度，不准确但是降低了计算量。

• 蒙特卡洛是一类随机算法，除此以外还有很多随机算法，比如拉斯维加斯算法（结果总是正确的算法）

## x. 参考教程

posted @ 2022-07-06 18:11  climerecho  阅读(1074)  评论(0编辑  收藏  举报