Методы ресэмплинга и Монте-Карло#

Введение#

Методы повторной выборки и Монте-Карло — это статистические техники, которые заменяют математический анализ большим количеством вычислений.

Например, представьте, что вы и ваш брат Кайл оказались автостопом на длинной и одинокой дороге. Внезапно появляется сияющий демон... посреди... дороги. И он говорит:

Если вы подбрасываете монету с вероятностью выпадения орла \(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.