您当前的位置:首页 > 电脑百科 > 程序开发 > 算法

蒙特卡洛算法

时间:2022-07-08 09:42:24  来源:博客园  作者:海椰人

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/ 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/

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()
折叠 

效果如图:

蒙特卡洛算法

 

以上内容参考:

  1. 课程视频
  2. https://www.section.io/engineering-education/buffon-needle/
  3. https://Github.com/topics/buffon-needle
  4. https://github.com/GunnarDahm/buffon_monte_carlo_sim/blob/master/buffon_monte_carlo.py
  5. https://blog.csdn.NET/qq_45757739/article/detAIls/108387567
  6. https://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



Tags:算法   点击:()  评论:()
声明:本站部分内容及图片来自互联网,转载是出于传递更多信息之目的,内容观点仅代表作者本人,不构成投资建议。投资者据此操作,风险自担。如有任何标注错误或版权侵犯请与我们联系,我们将及时更正、删除。
▌相关推荐
诱导付费、自动扣费……微短剧被质疑借助算法精准“围猎”老年人
诱导付费、自动扣费、重复收费&hellip;&hellip;聚焦身边的消费烦心事⑦丨一些微短剧被质疑借助算法精准“围猎”老年人中工网北京3月31日电(工人日报&mdash;中工网记者刘兵)...【详细内容】
2024-04-01  Search: 算法  点击:(11)  评论:(0)  加入收藏
分析网站SEO快速排名算法对网站具体的影响效果
亲爱的朋友们,今天我想和大家分享一个我们都关心的话题&mdash;&mdash;网站SEO快速排名算法对网站我们身处一个信息爆炸的时代,如何在海量的信息中脱颖而出,成为了一个我们不得...【详细内容】
2024-03-28  Search: 算法  点击:(21)  评论:(0)  加入收藏
当prompt策略遇上分治算法,南加大、微软让大模型炼成「火眼金睛」
近年来,大语言模型(LLMs)由于其通用的问题处理能力而引起了大量的关注。现有研究表明,适当的提示设计(prompt enginerring),例如思维链(Chain-of-Thoughts),可以解锁 LLM 在不同领域的...【详细内容】
2024-03-12  Search: 算法  点击:(21)  评论:(0)  加入收藏
谷歌宣布更新搜索算法:打击AI生成内容,提高搜索结果质量
IT之家 3 月 6 日消息,谷歌于当地时间 5 日发文宣布,针对用户对搜索结果质量下降的反馈,将对算法进行调整,旨在打击 AI 生成的内容以及内容农场等垃圾信息,使用户能够看到更多“...【详细内容】
2024-03-06  Search: 算法  点击:(44)  评论:(0)  加入收藏
小红书、视频号、抖音流量算法解析,干货满满,值得一看!
咱们中国现在可不是一般的牛!网上的网友已经破了十个亿啦!到了这个互联网的新时代,谁有更多的人流量,谁就能赢得更多的掌声哦~抖音、小红书、、视频号,是很多品牌必争的流量洼地...【详细内容】
2024-02-23  Search: 算法  点击:(18)  评论:(0)  加入收藏
雪花算法详解与Java实现:分布式唯一ID生成原理
SnowFlake 算法,是 Twitter 开源的分布式 ID 生成算法。其核心思想就是:使用一个 64 bit 的 long 型的数字作为全局唯一 ID。在分布式系统中的应用十分广泛,且 ID 引入了时间戳...【详细内容】
2024-02-03  Search: 算法  点击:(54)  评论:(0)  加入收藏
简易百科之什么是搜索引擎的PageRank算法?
简易百科之什么是搜索引擎的PageRank算法?在互联网时代,搜索引擎是我们获取信息的重要工具。而PageRank算法则是搜索引擎的核心技术之一,它决定了网页在搜索结果中的排名。那么...【详细内容】
2024-01-24  Search: 算法  点击:(57)  评论:(0)  加入收藏
PageRank算法揭秘:搜索引擎背后的魔法师的工作原理
PageRank(PR)算法是由谷歌创始人之一的拉里&middot;佩奇LarryPage命名的一种衡量网站页面重要性的方法。根据谷歌的说法,PageRank通过计算页面链接的数量和质量来粗略估计分...【详细内容】
2024-01-23  Search: 算法  点击:(46)  评论:(0)  加入收藏
程序开发中常用的十种算法,你用过几种?
当编写程序时,了解和使用不同的算法对解决问题至关重要。以下是C#中常用的10种算法,每个算法都伴随着示例代码和详细说明。1. 冒泡排序 (Bubble Sort):冒泡排序是一种简单的比...【详细内容】
2024-01-17  Search: 算法  点击:(46)  评论:(0)  加入收藏
百度最新的搜索引擎算法是什么样的?
百度搜索引擎算法是百度用来决定网页排名的算法。它是百度搜索技术的核心,也是百度作为全球最大的中文搜索引擎的基石。随着互联网的发展和用户需求的不断变化,百度搜索引擎算...【详细内容】
2024-01-10  Search: 算法  点击:(92)  评论:(0)  加入收藏
▌简易百科推荐
小红书、视频号、抖音流量算法解析,干货满满,值得一看!
咱们中国现在可不是一般的牛!网上的网友已经破了十个亿啦!到了这个互联网的新时代,谁有更多的人流量,谁就能赢得更多的掌声哦~抖音、小红书、、视频号,是很多品牌必争的流量洼地...【详细内容】
2024-02-23  二手车小胖说    Tags:流量算法   点击:(18)  评论:(0)  加入收藏
雪花算法详解与Java实现:分布式唯一ID生成原理
SnowFlake 算法,是 Twitter 开源的分布式 ID 生成算法。其核心思想就是:使用一个 64 bit 的 long 型的数字作为全局唯一 ID。在分布式系统中的应用十分广泛,且 ID 引入了时间戳...【详细内容】
2024-02-03   一安未来  微信公众号  Tags:雪花算法   点击:(54)  评论:(0)  加入收藏
程序开发中常用的十种算法,你用过几种?
当编写程序时,了解和使用不同的算法对解决问题至关重要。以下是C#中常用的10种算法,每个算法都伴随着示例代码和详细说明。1. 冒泡排序 (Bubble Sort):冒泡排序是一种简单的比...【详细内容】
2024-01-17  架构师老卢  今日头条  Tags:算法   点击:(46)  评论:(0)  加入收藏
百度推荐排序技术的思考与实践
本文将分享百度在推荐排序方面的思考与实践。在整个工业界的推广搜场景上,特征设计通常都是采用离散化的设计,需要保证两方面的效果,一方面是记忆,另一方面是泛化。特征都是通过...【详细内容】
2024-01-09  DataFunTalk  微信公众号  Tags:百度推荐   点击:(81)  评论:(0)  加入收藏
什么是布隆过滤器?如何实现布隆过滤器?
以下我们介绍了什么是布隆过滤器?它的使用场景和执行流程,以及在 Redis 中它的使用,那么问题来了,在日常开发中,也就是在 Java 开发中,我们又将如何操作布隆过滤器呢?布隆过滤器(Blo...【详细内容】
2024-01-05  Java中文社群  微信公众号  Tags:布隆过滤器   点击:(94)  评论:(0)  加入收藏
面向推荐系统的深度强化学习算法研究与应用
随着互联网的快速发展,推荐系统在各个领域中扮演着重要的角色。传统的推荐算法在面对大规模、复杂的数据时存在一定的局限性。为了解决这一问题,深度强化学习算法应运而生。本...【详细内容】
2024-01-04  数码小风向    Tags:算法   点击:(106)  评论:(0)  加入收藏
非负矩阵分解算法:从非负数据中提取主题、特征等信息
非负矩阵分解算法(Non-negativeMatrixFactorization,简称NMF)是一种常用的数据分析和特征提取方法,主要用于从非负数据中提取主题、特征等有意义的信息。本文将介绍非负矩阵分解...【详细内容】
2024-01-02  毛晓峰    Tags:算法   点击:(75)  评论:(0)  加入收藏
再谈前端算法,你这回明白了吗?
楔子 -- 青蛙跳台阶一只青蛙一次可以跳上一级台阶,也可以跳上二级台阶,求该青蛙跳上一个n级的台阶总共需要多少种跳法。分析: 当n=1的时候,①只需要跳一次即可;只有一种跳法,即f(...【详细内容】
2023-12-28  前端爱好者  微信公众号  Tags:前端算法   点击:(114)  评论:(0)  加入收藏
三分钟学习二分查找
二分查找是一种在有序数组中查找元素的算法,通过不断将搜索区域分成两半来实现。你可能在日常生活中已经不知不觉地使用了大脑里的二分查找。最常见的例子是在字典中查找一个...【详细内容】
2023-12-22  小技术君  微信公众号  Tags:二分查找   点击:(81)  评论:(0)  加入收藏
强化学习算法在资源调度与优化中的应用
随着云计算和大数据技术的快速发展,资源调度与优化成为了现代计算系统中的重要问题。传统的资源调度算法往往基于静态规则或启发式方法,无法适应动态变化的环境和复杂的任务需...【详细内容】
2023-12-14  职场小达人欢晓    Tags:算法   点击:(169)  评论:(0)  加入收藏
站内最新
站内热门
站内头条