6. 蒙特卡洛算法
6.1 计算π" role="presentation" style="display: inline; font-style: normal; font-weight: normal; text-indent: 0px; text-align: left; text-transform: none; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; white-space: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">ππ
回到顶部
a. 原理
如果我们不知道 ππ 的值，我们能不能用随机数 来近似 ππ 呢？
假设我们用一个随机数生成器，每次生成两个范围在 [−1,+1][−1,+1] 的随机数，一个作为 x，另一个作为 y，即生成了一个二维随机点：
假如生成 1亿 个随机样本，会有多少落在 半径=1 的圆内？这个概率就是圆的面积除以正方形的面积。
即：P=πr222=π4P=πr222=π4
假设从正方形区域中随机抽样 n 个点，那么落在圆内点个数的期望为：Pn=πn4Pn=πn4，
下面我们去求落在圆内的点的个数，只需满足x2+y2⩽1x2+y2⩽1 即为圆内。
如果生成的随机点的个数足够多，落在圆内的实际观测值 m≈πn4m≈πn4；
我们已知了m 与 n，所以π≈4mnπ≈4mn.
事实上，根据概率论大数定律：
4mn→π4mn→π，as n → ∞
这保证了蒙特卡洛的正确性。
伯恩斯坦概率不等式还能确定 观测值和真实值之间误差的上界。
|4mn−π|=O(1n√)|4mn−π|=O(1n)
说明 这个误差与样本n的根号成反比。
回到顶部
b. 代码
下面放一个Python代码
hide code#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
另外可还有一个可视化程序，可以模拟点落在方块区域圆内外：http://www.anders.wang/monte-carlo/
6.2 Buffon's Needle Problem
回到顶部
a. 原理
布封投针，也是用蒙特卡洛来近似 ππ 值。这是一个可以动手做的实验。
用一张纸，画若干等距平行线（距离为 d），撒上一把等长的针（长度为l），通过与平行线相交的针的数量，就可以推算出 ππ。
通过微积分可以算出：相交的概率为：P=2lπdP=2lπd
微积分推导过程：
课程里并没有讲解推导，这里我参考的是一下两篇博客的推导过程：
https://zhuanlan.zhihu.com/p/479953215https://cosx.org/2009/11/a-brief-talk-on-buffon-throwing-needle-problems/
主流做法是通过对针的斜率进行积分：
这里我后续补充。
跟 6.1 类似，我们随机扔 n 根针，这样相交个数的期望为 Pn=2lnπdPn=2lnπd 。我们可以观察到（如果是电脑模拟即为通过公式判断出）有 m 跟针实际与线相交，如果n足够大，则 m≈2lnπdm≈2lnπd。
求 ππ 公式即为： π≈2lnmdπ≈2lnmd
回到顶部
b. 代码
有了公式 π≈2lnmdπ≈2lnmd，代码实现其实很简单了，仅列出一种实现思路：
hide codeimport 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
当然，也有可视化的代码：
hide codeimport 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)]) self.end_points = np.array([np.add(self.needle_coordinates, -1*np.array( self.complex_representation)), np.add(self.needle_coordinates, self.complex_representation)]) 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() 折叠
效果如图：
以上内容参考：
课程视频https://www.section.io/engineering-education/buffon-needle/https://github.com/topics/buffon-needlehttps://github.com/GunnarDahm/buffon_monte_carlo_sim/blob/master/buffon_monte_carlo.pyhttps://blog.csdn.net/qq_45757739/article/details/108387567https://blog.csdn.net/TSzero/article/details/111604960理解思想即可，如果后续有机会，可能单出一篇介绍介绍，也有可能将这部分丰富一下。
6.3 估计阴影部分的面积
我们稍微推广一下，试着用蒙特卡洛解决一个阴影部分面积的求解。比如下图：
我们如何使用蒙特卡洛的思路解决这个阴影部分面积的求解呢？
类似于上面的思路，在正方形内做随机均匀抽样，得到很多点，怎么确定点在阴影部分呢？
可知，阴影部分的点满足：
{x2+y2>4(x−1)2+(y−1)2≤1{x2+y2>4(x−1)2+(y−1)2≤1
易知，正方形面积 A1=4A1=4；设阴影部分面积为 A2A2随机抽样的点落在阴影部分的概率为：P=A2A1=A24P=A2A1=A24从正方形区域抽样 n 个点，n尽可能大，则来自阴影部分点的期望为：nP=nA24nP=nA24；如果实际上满足上述条件的点 有 m 个，则令 m≈nPm≈nP得到：A2≈4mnA2≈4mn代码与 6.1 相近。
6.4 求不规则积分
近似求积分是蒙特卡洛在工程和科学问题中最重要的应用。很多积分是没有解析的积分（即可以计算出来的积分），特别是多元积分，而只能用数值 *** 求一个近似值，蒙特卡洛就是最常用的数值 *** 。
一元函数步骤如下：
我们要计算一个一元函数的定积分 I=∫baf(x)dxI=∫abf(x)dx;
从区间 [a,b][a,b] 上随机均匀抽样 x1,x2,...,xnx1,x2,...,xn;计算 Qn=(b−a)1n∑ni=1f(xi)Qn=(b−a)1n∑i=1nf(xi)，即均值乘以区间长度；这里均值乘以区间长度是 实际值，而 I 是期望值用 QnQn 近似 II大数定律保证了 当n→∞,Qn→In→∞,Qn→I
多元函数步骤如下：
我们要计算一个多元函数的定积分 I=∫baf(x⃗ )dx⃗ I=∫abf(x→)dx→，积分区域为 ΩΩ;
从区间 ΩΩ 上随机均匀抽样 x1→,x2→,...,xn→x1→,x2→,...,xn→;计算 ΩΩ 的体积V（高于三维同样）：V=∫Ωdx⃗ V=∫Ωdx→；hh值得注意的是，这一步仍要计算定积分，如果形状过于复杂，无法求得 V，那么无法继续进行，则无法使用蒙特卡洛算法。所以只能适用于比较规则的区域，比如圆形，长方体等。计算 Qn=V1n∑ni=1f(xi→)Qn=V1n∑i=1nf(xi→)，即均值乘以区间长度；这里均值乘以区间长度是 实际值，而 I 是期望值用 QnQn 近似 II下面我们从积分的角度再来看看 蒙特卡洛近似求 pi
定义一个二元函数 f(x,y)={1 if点在圆内0 if点在圆外f(x,y)={1 if点在圆内0 if点在圆外;定义一个区间 Ω=[−1,1]×[−1,1]Ω=[−1,1]×[−1,1]I=πr2=πI=πr2=π接下来用蒙特卡洛近似 I，得到关于 ππ的算式即可得到近似的ππ;随机抽样 n 个点，记为(x1,y1),...,(xn,yn)(x1,y1),...,(xn,yn)计算 区域面积 V=∫Ωdxdy=4V=∫Ωdxdy=4；计算 Qn=V1n∑ni=1f(xi,yi)Qn=V1n∑i=1nf(xi,yi)蒙特卡洛近似 Q 与 I 近似相等：π=Qn=∫Ωf(x,y)dxdyπ=Qn=∫Ωf(x,y)dxdy这是从蒙特卡洛积分的角度得到的pi，6.1 中则是从蒙特卡洛概率和期望的角度得到的。
6.5 用蒙特卡洛近似期望
这个 *** 对于统计学和机器学习很有用。
定义 X 是 d 维的随机变量，函数 p(x) 是一个PDF，概率密度函数；函数 f(x)f(x) 的期望：Ex∼p[f(X)]=∫Rdf(X)⋅(x)dxEx∼p[f(X)]=∫Rdf(X)⋅(x)dx直接以上面的方式求期望可能并不容易，所以通常使用蒙特卡洛近似求期望：随机抽样：根据概率密度函数 p(x)p(x) 进行随机抽样，记为X1,X2,...,XnX1,X2,...,Xn；计算 Qn=1n∑ni=1f(xi)Qn=1n∑i=1nf(xi)用 Q 近似 期望Ex∼p[f(X)]Ex∼p[f(X)]6.6 总结 | 蒙特卡洛算法的思想
我的想法是尽量精简，即：
模拟---抽样---估值，通过模拟出来的大量样本集或者随机过程，以随机抽样的方式，去近似我们想要研究的实际问题对象。
补充蒙特卡洛相关：
蒙特卡洛是摩洛哥的赌场；
蒙特卡洛算法得到的结果通常是错误的，但很接近真实值，对于对精度要求不高的机器学习已经足够。随机梯度下降就是一种蒙特卡洛算法，用随机的梯度近似真实的梯度，不准确但是降低了计算量。
蒙特卡洛是一类随机算法，除此以外还有很多随机算法，比如拉斯维加斯算法（结果总是正确的算法）
文章来自https://www.cnblogs.com/Roboduster/p/16451932.html
