重采样与蒙特卡洛方法#
简介#
重采样和蒙特卡洛方法是统计技术,它们用大量的计算代替数学分析。
例如,假设你和你的兄弟凯尔正在一条漫长而孤独的路上搭便车。突然,一个闪亮的恶魔出现在……路中间。他说:
如果你抛掷一枚正面朝上概率为 \(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
中的重采样和蒙特卡洛方法。