布封投针问题（Buffon's needle problem）

蒙特卡洛方法（Monte-Carlo Simulation）

估算PI

估计不规则图形的面积

随机抛点

采样估计

样本采集

拒绝采样（reject sample）

布封投针问题（Buffon's needle problem）

问题：

1、取一张白纸，画出间隔为 $d$ 的等距平行线。

2、取一根长度为 $l(l\leq d)$ 的针，随机地向画有等距平行线的纸上掷 $n$ 次，观察针与直线相交的次数，记为 $m$ 。

3、计算针与直线相交的概率。

这个概率值是 $p=\frac{2l}{\pi d}$ ，其中 $\pi$ 是圆周率。

证明：

上图中， $x$ 表示针的尾部与平行线的距离， $x\in [0,d]$ ， $\theta$ 表示针与平行线的夹角， $\theta \in [0,\pi ]$ 。由于针是随机落到平行线内的，所以概率是均等的，即 $\frac{1}{\pi d}$ 。要使得针与边界相交，那么就要满足 $x\leq lsin\theta$ 。那么要求针与边界相交的概率，就是在可行域内解进行积分，即：

$p=\int_{0}^{\pi }\int_{0}^{lsin\theta }\frac{1}{\pi d}dxd\theta \Rightarrow \frac{2l}{\pi d}$

Buffon投针问题，反映出了，相交次数与 $\pi$ 的关系，从而可以估算出 $\pi$ ，Buffon投针是首次从几何角度来表达概率的，也是首次通过随机试验来处理一个确定性的数学问题。

蒙特卡洛方法（Monte-Carlo Simulation）

蒙特卡洛方法是一种以概率统计理论为指导的一类非常重要的数值计算方法。蒙特卡罗是摩纳哥公国的一座城市（赌城）。Buffon投针就是蒙特卡洛方法的思想，但是Buffon是蒙特卡洛方法的起源。蒙特卡洛方法同样可以通过随机通过产生随机数的方式来解决计算问题。比如估算 $\pi$ 和计算不规则图形的面积。如果你觉得Buffon投针方法来估算 $\pi$ 比较复杂，那么来看用另一种简单的Monte-Carlo方法估算 $\pi$ 。

估算PI

给定一个边长为 $1$ 的正方形，里是一个内切的 $\frac{1}{4}$ 圆，我们知道正方形的面积是 $1$ ，而 $\frac{1}{4}$ 的面积 $\frac{\pi }{4}$ 。怎么估算 $\pi$ 呢？向正方形内随机抛点，最后计算落在圆弧内的点，计算概率然后再乘以 $4$ ，就是 $\pi$ 的估算值。如下图所示。随着抛出的点数增加， $\pi$ 的估计值越准确。

# -*- coding: utf-8 -*-
"""
Created on Sat Sep  1 11:07:51 2018@author: abner_hg
"""import random
import mathicount = 0
iter_num = 10000
x_axis = []
y_axis = []
for i in xrange(iter_num):x = random.random()y = random.random()if math.sqrt(x**2 + y**2) <= 1:icount += 1
print('Pi = ', float(4*icount)/float(iter_num))

估计不规则图形的面积

随机抛点估计

假设我们需要计算一个不规则图形的面积，那么图形的不规则程度和分析性计算（比如积分）的复杂程度是成正比的。而采用蒙特卡罗方法是怎么计算的呢？如下图，假设我们不知道曲线的具体函数表达式或者难以求出具体的函数表达式，那么我们该怎么估算阴影部分的面积？首先可以取一个边长为 $a$ 的正方形，然后你可以随机地向正方形内，随机抛一些点，计算落在阴影部分的个数，最后除以抛出的总个数，得出的值再乘以 $a^{2}$ ，就是阴影部分的面积一个估算。抛出的点数越多，估算越准确，这个就是蒙特卡洛思想。

采样估计

假设我们知道函数表达式 $f(x)$ ，我们知道求解阴影部分的面积，可以对函数直接进行积分 $\int_{0}^{a}f(x)dx$ 。那么如果，我们难以计算积分，该怎么办？随机抛点的方法是可以的，但是有没有其他的方法？如下图，我们可以随机地在 $[0,a]$ 上取一个点 $x_{0}$ ，用 $f(x_{0})$ 来代表所有的函数值，那么积分值就可以认为是 $(a-0)f(x_{0})=af(x_{0})$ 。

上图的估算结果肯定是不准确的，因为从图中来看，有很多函数值是大于 $f(x_{0})$ 的，那么为了使得估算更精确，我们需要采集更多的点，如下图，

以这些函数值的均值 $\frac{f(x_{0})+\cdots +f(x_{n})}{n}$ 来表示所有的函数值,那么对积分的估算，就变成了 $\frac{(a-0)[f(x_{0})+\cdots +f(x_{n})]}{n}$ 。这种方式估算貌似比上一种更靠谱些。但是存在一个问题，我们取函数值都是随机取的，也就是在哪个地方采集样本点的概率都是一样的。假如我们的函数 $f(x_{0})$ 表示概率密度函数（probability density function，pdf），积分就变成了求其累积分布函数（Cumulative Distribution Function，CDF）。那么在上图中 $x_{i}$ 附近的点，应该比 $x_{k}$ 附近的点更加密集，因为 $x_{i}$ 附近处的概率值比较大，而 $x_{k}$ 附近处的概率值比较小，因此我们不能简单地随机采集样本。那么求解上述积分的关键，就是如何采集样本点。

样本采集

如何采集样本是蒙特卡洛方法的关键一步，在《LDA数学八卦》中，有这么一段话：

统计模拟中有一个重要的问题就是给定一个概率分布 $p(x)$ ，我们如何在计算机中生成它的样本。一般而言均匀分布 $Uniform(0,1)$ 的样本是相对容易生成的。通过线性同余发生器可以生成伪随机数，我们用确定性算法生成[0,1]之间的伪随机数后，这些序列的各种统计指标和均匀分布 $Uniform(0,1)$ 的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。

我们常见的概率分布，连续或者离散的情况，都可以基于 $Uniform(0,1)$ 的样本生成。例如正态分布，可以由Box-Muller变换得到。

定理：（Box-Muller变换）如果随机变量 $U_{1}$ ， $U_{2}$ 独立且 $U_{1},U_{2}\sim Uniform[0,1]$ ，则：

$Z_{0}=\sqrt{-2lnU_{1}}cos(2\pi U_{2})$

$Z_{1}=\sqrt{-2lnU_{1}}sin(2\pi U_{2})$

则 $Z_{1},Z_{2}$ 独立且服从标准正态分布。

写这么多的目的就是想表达，很多常见的分布，是可以基于均匀分布来实现采样。回到刚才的问题，我们的 $f(x)$ 比较复杂，难以采样，那么该如何实现采样？

拒绝采样（reject sample）

接下来我们以概率分布函数来讲解，假如给定的函数的值域大于1，我们可以归一化到 $[0,1]$ 之间，从而可以看成概率密度函数。假如 $q(x)$ 是我们常见的可采集的概率密度函数，那么是否存在一个 $q(x)$ 永远大于 $f(x)$ 呢？答案是否定的，因为所有的概率之和要满足1，故任意两个概率密度函数一定有交点。如下图，就是我们找到的一个概率密度函数。

拒绝采样的步骤：

1、在 $q(x)$ 上采集一个样本点 $x_{i}$

2、计算接受率 $\alpha=\frac{f(x_{i})}{mq(x_{i})}$ 。

3、在均匀分布 $u\sim Uniform(0,1)$ 上，随机选择一个 $u$ 。

4、如果 $u\leq \alpha$ 那么接受 $x_{i}$ 。

5、重复步骤 $1\sim 4$ 。

解释两个问题，1、什么是接受率。2、为什么要在均匀分布上选取一个 $u$ 来决定是否接受 $x_{i}$ 。

问题1：

如上图所示，表示落在绿线以下的部分是可接受的，落在绿线和红线之间的，是被拒绝的。那么那个判断的临界点就是接受率，在数学上的表示就是 $\alpha=\frac{f(x_{i})}{mq(x_{i})},\alpha \in [0,1]$ 。

问题2：

我们在 $q(x)$ 上随机采集了一个点 $x_{0}$ ，我们不知道它是在那个区域，我们凭直觉可以知道，如果两个函数，越接近，那么被接受的区域就越大（上图紫线越长），也就是 $\alpha$ 越大，如果两个函数相聚越远，那么被拒绝的区域也就越大（上图红线越长），也就是 $\alpha$ 越小。此时怎么办，不能而用抛硬币的方法来决断，太过草率，由于 $\alpha$ 是一个属于 $[0,1]$ 之间的一个值，我们可以用在均匀分布上进行随机生成一个数，如果大于 $\alpha$ 就拒绝，如果小于就接受。当然，这种方式也存在误判情况。但是如果你的 $q(x)$ 与 $f(x)$ 比较接近的话，被误判的可能性就越小，但是随着采集的样本点的增加，被接受的样本点的函数值，大体上会与 $f(x)$ 保持一致。

拒绝采样的缺点：

1、如果 $mq(x)$ 与 $f(x)$ 相距比较远，那么被被拒绝的样本点会很多，采集10000个有可能被拒绝掉8000个，即使相距较劲，也可能被拒绝掉，费了很大劲，才计算出来的结果，被拒绝了，费时费力！

2、如果我们的分布是在高维上， $q(x)$ 与 $m$ 是很难确定的。

采样的方法有很多，以后有时间会在这一块写个总结。

我们知道了如何采样，就知道了样本集，那么我们就可以对计算问题利用随机模拟的方式进行估算，这就是蒙特卡洛方法的思想。

参考：《LDA数学八卦》

蒙特卡洛方法（Monte-Carlo Simulation）相关推荐

Stata: 蒙特卡洛模拟(Monte Carlo Simulation)没那么神秘
作者:侯新烁湘潭大学 [编译] (知乎 | 简书 | 码云) Stata连享会精彩推文1 || 精彩推文2 资料参考来源: The Stata Blog » Monte Carlo simulat ...
心得复述知识体系：《强化学习》中的蒙特卡洛方法 Monte Carlo Methods in Reinforcement Learning
前言: 刚刚读完 Sutton 的<强化学习(第二版)>第5章:蒙特卡洛方法.为了巩固本章收获,笔者将在本文中用尽量简单直白的语言复述本章的思想,各个知识点之间的关系.同时,这方便笔者日后 ...
【路径追踪】数学工具--蒙特卡洛方法(Monte Carlo)
Intro 蒙特卡洛方法是一类通过随机采样来求解问题的算法, 要求解的问题是某随机事件的概率或某随机变量的期望. 现在认为最早记载的一个蒙特卡洛计算示例是由蒲丰在 1777 年完成的投针试验. 在实验 ...
java计算椭圆的面积_java算法3_蒙特卡洛方法(Monte Carlo method)求PI和椭圆面积
蒙特卡洛方法,是一种以概率统计理论为指导的一类非常重要的数值计算方法.是指使用随机数来解决很多计算问题的方法.蒙特卡洛方法的名字来源于摩纳哥的一个城市蒙特卡洛,该城市以×××业闻名,而蒙特卡洛方法正是 ...
【Matlab 六自由度机器人】基于蒙特卡罗方法(Monte Carlo Method)构建机器人工作空间（附MATLAB建模仿真完整代码）
[Matlab 六自由度机器人]基于蒙特卡洛法构建工作空间往期回顾前言正文一.蒙特卡洛方法及机器人工作空间的概述 1. 蒙特卡洛法 2. 机器人工作空间二.rand函数的应用 1. rand ...
correlation 蒙特卡洛_蒙特卡洛模拟法及其matlab案例（Monte Carlo simulation method and its matlab case）.doc...
蒙特卡洛模拟法及其matlab案例(Monte Carlo simulation method and its matlab case) 蒙特卡洛模拟法及其matlab案例(Monte Carlo s ...
蒙特卡洛（Monte Carlo）法求定积分
蒙特卡洛(Monte Carlo)法是一类随机算法的统称.随着二十世纪电子计算机的出现,蒙特卡洛法已经在诸多领域展现出了超强的能力.在机器学习和自然语言处理技术中,常常被用到的MCMC也是由此发展而来 ...
蒙特卡洛积分(Monte Carlo Integration)应用：利用蒙特卡洛积分生成 McBeth表
蒙特卡洛积分(Monte Carlo Integration)应用蒙特卡洛积分通常函数f(x)的积分: 可以解释为计算函数曲线下方的面积: 而我们的蒙特卡洛积分则是通过近似的方式来获取一个函数的积 ...
matlab monte carlo,Monte Carlo Simulation
Monte Carlo simulation is a technique used to study how a model responds to randomly generated input ...
MATLAB蒙特卡洛（Monte Carlo）方法求椭圆面积
MATLAB蒙特卡洛方法求椭圆面积代码代码在某个规定的范围内随机打点,找到满足条件的点,并数一下这些点的数量与总的随机点数量的比,就OK了.关键是设置条件. 代码 clear;clc; n=10 ...

蒙特卡洛方法（Monte-Carlo Simulation）

布封投针问题（Buffon's needle problem）

蒙特卡洛方法（Monte-Carlo Simulation）

估算PI

估计不规则图形的面积

随机抛点估计

采样估计

样本采集

拒绝采样（reject sample）

蒙特卡洛方法（Monte-Carlo Simulation）相关推荐

最新文章

热门文章