回到顶部
如果我们不知道 ππ 的值,我们能不能用随机数 来近似 ππ 呢?
假设我们用一个随机数生成器,每次生成两个范围在 [−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的根号成反比。
回到顶部
下面放一个Python/ target=_blank class=infotextkey>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/
回到顶部
布封投针,也是用蒙特卡洛来近似 ππ 值。这是一个可以动手做的实验。
用一张纸,画若干等距平行线(距离为 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
回到顶部
有了公式 π≈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()
折叠
效果如图:
以上内容参考:
理解思想即可,如果后续有机会,可能单出一篇介绍介绍,也有可能将这部分丰富一下。
我们稍微推广一下,试着用蒙特卡洛解决一个阴影部分面积的求解。比如下图:
我们如何使用蒙特卡洛的思路解决这个阴影部分面积的求解呢?
类似于上面的思路,在正方形内做随机均匀抽样,得到很多点,怎么确定点在阴影部分呢?
可知,阴影部分的点满足:
{x2+y2>4(x−1)2+(y−1)2≤1{x2+y2>4(x−1)2+(y−1)2≤1
代码与 6.1 相近。
近似求积分是蒙特卡洛在工程和科学问题中最重要的应用。很多积分是没有解析的积分(即可以计算出来的积分),特别是多元积分,而只能用数值方法求一个近似值,蒙特卡洛就是最常用的数值方法。
一元函数步骤如下:
我们要计算一个一元函数的定积分 I=∫baf(x)dxI=∫abf(x)dx;
大数定律保证了 当n→∞,Qn→In→∞,Qn→I
多元函数步骤如下:
我们要计算一个多元函数的定积分 I=∫baf(x⃗ )dx⃗ I=∫abf(x→)dx→,积分区域为 ΩΩ;
下面我们从积分的角度再来看看 蒙特卡洛近似求 pi
这是从蒙特卡洛积分的角度得到的pi,6.1 中则是从蒙特卡洛概率和期望的角度得到的。
这个方法对于统计学和机器学习很有用。
我的想法是尽量精简,即:
模拟---抽样---估值,通过模拟出来的大量样本集或者随机过程,以随机抽样的方式,去近似我们想要研究的实际问题对象。
补充蒙特卡洛相关:
蒙特卡洛是摩洛哥的赌场;
蒙特卡洛算法得到的结果通常是错误的,但很接近真实值,对于对精度要求不高的机器学习已经足够。随机梯度下降就是一种蒙特卡洛算法,用随机的梯度近似真实的梯度,不准确但是降低了计算量。
蒙特卡洛是一类随机算法,除此以外还有很多随机算法,比如拉斯维加斯算法(结果总是正确的算法)
文章来自
https://www.cnblogs.com/Roboduster/p/16451932.html