重采样和蒙特卡洛方法#

简介#

重采样和蒙特卡洛方法是统计技术，它们用大量的计算来代替数学分析。

例如，假设你和你的兄弟凯尔发现自己搭便车沿着一条漫长而孤独的道路行驶。突然，一个闪亮的恶魔出现在道路中间。他说：

如果你抛掷一枚正面概率为 \(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 中的重采样和蒙特卡洛方法。