重新取樣與蒙地卡羅方法#

簡介#

重新取樣與蒙地卡羅方法是統計技術，以大量計算取代數學分析。

舉例來說，假設你和你的兄弟凱爾發現自己搭便車走在一條又長又寂寞的道路上。突然，一個閃亮的惡魔…在路中央…發光。他說

如果你擲一枚正面機率為 \(p=0.5\) 的硬幣，正好 \(n=100\) 次，正面次數小於或等於 \(x=45\) 的機率是多少？正確回答，否則我就吃掉你們的靈魂。
>>> import math
>>> import numpy as np
>>> p = 0.5  # probability of flipping heads each flip
>>> n = 100  # number of coin flips per trial
>>> x = 45  # we want to know the probability that the number of heads per trial will be less than or equal to this

你的兄弟凱爾是個分析型的人。他回答

隨著擲硬幣次數增加，正面次數的分布將趨於常態分布，平均數為 \(\mu = p n\)，標準差為 \(\sigma = \sqrt{n p (1 - p)}\)，其中 \(p = 0.5\) 是正面的機率，\(n=100\) 是擲硬幣的次數。\(x=45\) 次正面的機率可以近似為此常態分布的累積密度函數 \(F(x)\)。具體來說

\[F(x; \mu, \sigma) = \frac{1}{2} \left[ 1 + \mbox{erf} \left( \frac{x-\mu}{\sigma \sqrt{2}} \right) \right]\]

>>> # Kyle's Analytical Approach
>>> mean = p*n
>>> std = math.sqrt(n*p*(1-p))
>>> # CDF of the normal distribution. (Unfortunately, Kyle forgets a continuity correction that would produce a more accurate answer.)
>>> prob = 0.5 * (1 + math.erf((x - mean) / (std * math.sqrt(2))))
>>> print(f"The normal approximation estimates the probability as {prob:.3f}")
The normal approximation estimates the probability as 0.159

你比較務實，所以你決定採取計算方法（或更精確地說，蒙地卡羅方法）：只需模擬許多擲硬幣的序列，計算每次投擲中正面的次數，並將機率估計為正面次數不超過 45 的序列比例。

>>> # Your Monte Carlo Approach
>>> N = 100000  # We'll do 100000 trials, each with 100 flips
>>> rng = np.random.default_rng()  # use the "new" Generator interface
>>> simulation = rng.random(size=(n, N)) < p  # False for tails, True for heads
>>> counts = np.sum(simulation, axis=0)  # count the number of heads each trial
>>> prob = np.sum(counts <= x) / N  # estimate the probability as the observed proportion of cases in which the count did not exceed 45
>>> print(f"The Monte Carlo approach estimates the probability as {prob:.3f}")
The Monte Carlo approach estimates the probability as 0.187

惡魔回答

你們兩個都不正確。這個機率由二項分布給出。具體來說。

\[\sum_{i=0}^{x} {n \choose i} p^i (1-p)^{n-i}\]

>>> # The Demon's Exact Probability
>>> from scipy.stats import binom
>>> prob = binom.cdf(x, n, p)
>>> print(f"The correct answer is approximately {prob}")
The correct answer is approximately 0.18410080866334788

當你的靈魂被吃掉時，你可以從你的簡單蒙地卡羅方法比常態近似更準確的知識中獲得安慰。這並不罕見：當確切答案未知時，計算近似通常比分析近似更準確。此外，惡魔很容易發明分析近似（更不用說確切答案）不可用的問題。在這種情況下，計算方法是唯一的出路。

重新取樣與蒙地卡羅方法教學#

雖然最好在可以使用確切方法時使用它，但學習使用計算統計技術可以提高依賴分析近似的 scipy.stats 功能的準確性，顯著擴展您的統計分析能力，甚至提高您對統計的理解。以下教學將幫助您開始使用 scipy.stats 中的重新取樣和蒙地卡羅方法。