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

hearthougan發表於2018-09-01

目錄

 

布封投針問題（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 math

icount = 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數學八卦》

相關文章