5 приемов использования Scipy.stats для моделирования сценариев «что если»

В этой статье мы заглянем под капот библиотеки scipy.stats и рассмотрим пять важных приемов для создания высокопроизводительных и точных симуляций с использованием только NumPy и SciPy.

# Введение

Данные редко бывают статичными. Решения редко бывают безрисковыми. Как специалисту по анализу данных, вам часто приходится проверять на прочность бизнес-предположения, исследовать неопределенность распределения или моделировать альтернативные реальности.

• «Что если наши ежедневные затраты на привлечение активных пользователей удвоятся?»
• «Что если во время рекламной акции трафик на нашем сервере резко возрастет на 300%?»
• «Какова вероятность того, что наши операционные убытки в этом квартале превысят 50 000 долларов?»

Для ответа на подобные вопросы типа «а что если» необходимо перейти от простых точечных оценок (например, простого среднего) к надежному вероятностному мышлению. Хотя многие специалисты сразу же могут обратиться к сложным системам моделирования, стандартный научный стек Python уже содержит недооцененный инструмент именно для такого рода моделирования: scipy.stats. Помимо своей распространенной репутации инструмента для проведения простых проверок гипотез или вычисления p-значений, scipy.stats предоставляет унифицированный программный интерфейс для параметризации, выборки и вычисления показателей риска для десятков непрерывных и дискретных вероятностных распределений.

В этой статье мы заглянем под капот библиотеки scipy.stats и рассмотрим пять важных приемов для создания высокопроизводительных и точных симуляций с использованием только NumPy и SciPy.

# 1. Замораживание распределений для параметризации сценариев

При моделировании сценариев часто возникает необходимость представить различные состояния мира: консервативный базовый сценарий, оптимистичный наилучший сценарий и пессимистичный наихудший сценарий. В стандартном процедурном коде это можно представить, используя словари параметров (например, location loc и scale scale) и распаковывая их в функции каждый раз, когда необходимо оценить вероятность или взять образец.

Превосходный объектно-ориентированный подход заключается в замораживании распределений. В библиотеке scipy.stats вызов класса распределения (например, stats.norm, stats.lognorm или stats.gamma) и прямая передача параметров конструктору возвращает «замороженную» случайную переменную (экземпляр rv_frozen).

Этот заблокированный объект инкапсулирует всю вероятностную модель. Вы можете передавать его по коду как единый объект, менять определения сценариев на лету и вызывать такие методы, как .rvs(), .pdf(), .cdf() или .ppf(), не указывая параметры повторно.

Предположим, мы моделируем ежедневный спрос на продукцию. При ручной реализации нам необходимо перетаскивать словари или переменные на каждом этапе нашего конвейера обработки:

import scipy.stats as stats # Определение исходных параметров сценария scenarios = { «recession»: {«mean»: 800, «std»: 250}, «baseline»: {«mean»: 1200, «std»: 150}, «boom»: {«mean»: 1800, «std»: 300} } # Для получения выборок или оценки риска требуется ручная распаковка параметров def simulate_demand(scenario_name, size=5): params = scenarios[scenario_name] # Передача параметров в каждый статистический вызов samples = stats.norm.rvs(loc=params[«mean»], scale=params[«std»], size=size) p_exceed_1000 = 1 — stats.norm.cdf(1000, loc=params[«mean»], scale=params[«std»]) return samples, p_exceed_1000 for name in scenarios.keys(): samples, prob = simulate_demand(name) print(f»{name.capitalize()} -> Prob > 1000: {prob:.2%}»)

Вот более элегантная альтернатива. Создав экземпляр дистрибутива, SciPy фиксирует параметры и предоставляет нам самодостаточный, чистый объект сценария:

import scipy.stats as stats # С помощью scipy параметры распределения замораживаются в именованный объект scenario_models = { «recession»: stats.norm(loc=800, scale=250), «baseline»: stats.norm(loc=1200, scale=150), «boom»: stats.norm(loc=1800, scale=300) } # Конвейер просто принимает замороженный объект случайной величины, отделяя логику от параметров def run_scenario_pipeline(rv_frozen, size=5): # Методы блокировки вызываются непосредственно на объекте samples = rv_frozen.rvs(size=size) # sf() — функция выживания (1 — CDF) p_exceed_1000 = rv_frozen.sf(1000) return samples, p_exceed_1000 for name, model in scenario_models.items(): _, prob = run_scenario_pipeline(model) print(f»{name.capitalize()} Scenario | Prob demand > 1000: {prob:.2%}»)

Выход:

Сценарий рецессии | Вероятность спроса > 1000: 21,19% Базовый сценарий | Вероятность спроса > 1000: 90,88% Сценарий подъема | Вероятность спроса > 1000: 99,62%

Замораживая параметры, мы отделяем математические предположения от логики выполнения. Если мы хотим переключить нашу модель спроса с нормального распределения на асимметричное логарифмически нормальное распределение, нам нужно изменить всего одну строку, где мы инициализируем замороженный объект. Дальнейший конвейер обработки остается неизменным.

# 2. Моделирование методом Монте-Карло с использованием функции .rvs() для количественной оценки неопределенности

В бизнес-планировании специалисты часто создают электронные таблицы, в которых доход рассчитывается на основе статических точечных оценок — например, доход = ежедневный трафик * коэффициент конверсии * средняя стоимость заказа. Однако каждый из этих параметров обладает высокой степенью неопределенности. Умножение средних значений скрывает накопительную дисперсию, что приводит к ошибке средних значений или систематической недооценке риска.

Для количественной оценки этой неопределенности мы можем использовать моделирование методом Монте-Карло. Вместо того чтобы предполагать статичные значения, мы представляем каждую переменную в виде распределения вероятностей. Используя метод .rvs(size=n) на наших замороженных распределениях, мы можем мгновенно получить $N$ случайных выборок из всех входных данных, обработать их с помощью нашей собственной формулы в векторизованном конвейере NumPy и оценить полное распределение вероятностей конечного результата.

Расчет статического наилучшего/наихудшего случая не учитывает совместную вероятность событий и фактическое распределение результатов, а написание циклов вручную на чистом Python невероятно медленно.

import random # Хрупкие точечные оценки avg_traffic = 50000 avg_conversion = 0.05 avg_aov = 60.0 expected_revenue = avg_traffic * avg_conversion * avg_aov print(f»Точечная оценка ожидаемой выручки: ${expected_revenue:,.2f}») # Медленный итеративный цикл ручной выборки n_sims = 100000 revenue_clunky = [] for _ in range(n_sims): # Моделирование нормального трафика и равномерных коэффициентов конверсии traffic = random.gauss(50000, 5000) conversion = random.uniform(0.04, 0.06) aov = random.gammavariate(15, 4.0) revenue_clunky.append(traffic * conversion * aov)

Выход:

Предварительная оценка ожидаемой выручки: 150 000,00 долларов США.

Используя модели распределения scipy.stats и одновременно отбирая выборки с помощью .rvs(), мы можем выполнить всю симуляцию за одну векторизованную операцию NumPy. Давайте решим эту задачу в 4 отдельных шага:

1. Определите подходящие формы распределения для нашего сценария.
2. Возьмите N образцов.
3. Распространение через бизнес-логику (посредством векторизации)
4. Извлечение процентилей риска

import numpy as np import scipy.stats as stats # 1. Определяем подходящие формы распределения для нашего сценария np.random.seed(42) n_simulations = 100_000 # Трафик симметричный (нормальный) traffic_dist = stats.norm(loc=50000, scale=5000) # Коэффициент конверсии — это вероятность, ограниченная значениями от 0 до 1 (бета) # Среднее значение = 5 / (5 + 95) = 5% conversion_dist = stats.beta(a=5, b=95) # Средняя стоимость заказа положительная и имеет правостороннюю асимметрию (гамма) # Среднее значение = 15 * 4 = 60 долларов aov_dist = stats.gamma(a=15, scale=4) # 2. Выбираем N выборок traffic_samples = traffic_dist.rvs(size=n_simulations) conversion_samples = conversion_dist.rvs(size=n_simulations) aov_samples = aov_dist.rvs(size=n_simulations) # 3. Векторизованное распространение через бизнес-логику revenue_samples = traffic_samples * conversion_samples * aov_samples # 4. Извлечение процентилей риска mean_rev = np.mean(revenue_samples) # 5% вероятность получить меньше p5_rev = np.percentile(revenue_samples, 5) # 5% вероятность получить больше p95_rev = np.percentile(revenue_samples, 95) print(f»Вероятностное среднее значение дохода: ${mean_rev:,.2f}») print(f»5-й процентиль (снижение): ${p5_rev:,.2f}») print(f»95-й процентиль (Вверху): ${p95_rev:,.2f}»)

Выход:

Вероятностное среднее значение выручки: 150 153,16 долл. США. 5-й процентиль (снижение): 51 294,37 долл. США. 95-й процентиль (рост): 300 734,30 долл. США.

Хотя средняя выручка соответствует нашей первоначальной точечной оценке (~150 тыс. долларов), моделирование методом Монте-Карло показывает, что существует 5%-ная вероятность того, что наша выручка упадет ниже 97 180 долларов из-за совокупной волатильности трафика, конверсии и стоимости заказа. Точечные оценки полностью скрывают этот операционный риск.

# 3. Анализ чувствительности с помощью параметрического анализа

В сценарном анализе вам может потребоваться выяснить, насколько чувствителен ваш риск снижения стоимости к изменениям конкретных исходных предположений. Например, если вы моделируете затраты на привлечение клиентов (CAC), вам нужно понять, как изменение волатильности маркетинга (стандартного отклонения) влияет на наихудший сценарий (95-й процентиль) CAC.

Хотя для каждой конфигурации можно было бы провести полномасштабное моделирование методом Монте-Карло, библиотека scipy.stats предлагает гораздо более быстрый и математически элегантный способ: функцию процентных точек (.ppf()), которая является обратной функцией кумулятивной функции распределения (CDF).

Путем изменения параметров, фиксации распределений и выполнения функции .ppf() мы можем аналитически рассчитать точные границы процентилей за микросекунды, не генерируя ни одной случайной выборки.

Моделирование тысяч выборок для каждой перестановки параметров только для того, чтобы найти процентиль, крайне неэффективно и требует больших вычислительных затрат.

import numpy as np import scipy.stats as stats mean_cac = 50.0 volatilities = [5.0, 10.0, 15.0, 20.0] # Запуск масштабной случайной симуляции только для извлечения одного процентиля for vol in volatilities: samples = stats.norm.rvs(loc=mean_cac, scale=vol, size=100000) p95_clunky = np.percentile(samples, 95) print(f»Std: {vol} -> 95th Percentile CAC: ${p95_clunky:.2f} (sampled)»)

Пример выходных данных:

Стандарт: 5,0 -> 95-й процентиль Стоимость привлечения клиента: 58,23 долл. США (по выборке) Стандарт: 10,0 -> 95-й процентиль Стоимость привлечения клиента: 66,53 долл. США (по выборке) Стандарт: 15,0 -> 95-й процентиль Стоимость привлечения клиента: 74,65 долл. США (по выборке) Стандарт: 20,0 -> 95-й процентиль Стоимость привлечения клиента: 82,85 долл. США (по выборке)

Используя метод .ppf() на наших замороженных распределениях, мы мгновенно вычисляем точный аналитический порог.

import scipy.stats as stats mean_cac = 50.0 volatilities = [5.0, 10.0, 15.0, 20.0] # Метод SciPy: перебор параметров и вычисление точных процентилей с помощью .ppf() for vol in volatilities: # 1. Замораживаем распределение cac_dist = stats.norm(loc=mean_cac, scale=vol) # 2. Вычисляем точный 95-й процентиль аналитически p95_exact = cac_dist.ppf(0.95) # 3. Вычисляем вероятность превышения экстремального порога в 80 долларов p_exceed_80 = cac_dist.sf(80.0) print(f»Волатильность (стандартное отклонение): ${vol:02.0f} | 95-й процентиль CAC: ${p95_exact:.2f} | P(CAC > $80): {p_exceed_80:.2%}»)

Выход:

Волатильность (стандартное отклонение): 0,05 долл. | 95-й процентиль CAC: 58,22 долл. | P(CAC > 80 долл.): 0,00% Волатильность (стандартное отклонение): 10 долл. | 95-й процентиль CAC: 66,45 долл. | P(CAC > 80 долл.): 0,13% Волатильность (стандартное отклонение): 15 долл. | 95-й процентиль CAC: 74,67 долл. | P(CAC > 80 долл.): 2,28% Волатильность (стандартное отклонение): 20 долл. | 95-й процентиль CAC: 82,90 долл. | P(CAC > 80 долл.): 6,68%

Этот анализ выявляет наши предельные возможности: если волатильность наших затрат на приобретение вырастет с 5 до 20 долларов, то 95-й процентиль стоимости приобретения подскочит с 58,22 до 82,90 долларов , а риск превышения нашего максимального бюджета на приобретение в 80 долларов резко возрастет с 0,00% до 6,68% .

# 4. Моделирование риска в хвостах распределения с помощью распределений с тяжелыми хвостами

Распространенная ошибка в сценарном планировании — предположение, что каждый непрерывный набор данных подчиняется стандартному гауссовому (нормальному) распределению. Хотя нормальное распределение легко моделировать, у него чрезвычайно узкие хвосты. В реальном мире простои системы, финансовые потери и задержки API имеют «тяжелые хвосты» : экстремальные события происходят гораздо чаще, чем может предсказать гауссовская модель.

Для надлежащего стресс-тестирования наших систем на предмет риска, связанного с «хвостовыми» распределениями, мы должны заменить наивные предположения о нормальном распределении альтернативами с «тяжелыми хвостами», такими как t-критерий Стьюдента (stats.t), логнормальное распределение (stats.lognorm) или распределение Парето (stats.pareto).

Используя метод `.fit()` из библиотеки `scipy.stats`, мы можем аппроксимировать исторические наблюдения как нормальным распределением, так и распределением с «тяжелыми хвостами», а затем использовать функцию выживания (`.sf()`) для количественной оценки реальной вероятности наихудшего сценария отказов.

При работе с данными, имеющими «тяжелые хвосты» распределения, моделирование с использованием нормального распределения полностью лишает вас возможности отслеживать экстремальные события с «нижними» хвостами:

import numpy as np import scipy.stats as stats # Синтетические исторические данные о задержке с тяжелыми хвостами (критерий Стьюдента с df=3) np.random.seed(42) latencies = stats.t(df=3, loc=200, scale=40).rvs(size=1000) # Наивно предполагая нормальное распределение без проверки соответствия mean_lat, std_lat = np.mean(latencies), np.std(latencies) prob_outage_clunky = 1 — stats.norm.cdf(400, loc=mean_lat, scale=std_lat) print(f»Наивное гауссовское P(Задержка > 400 мс): {prob_outage_clunky:.6%}»)

Выход:

Наивный гауссовский P (задержка > 400 мс): 0,046321%

Путем сопоставления распределения Стьюдента с нормальным распределением мы можем явно сравнить качество подгонки и точно рассчитать риски в хвостах распределения.

import numpy as np import scipy.stats as stats # Синтетические исторические данные о задержке (с тяжелым хвостом) np.random.seed(42) latencies = stats.t(df=3, loc=200, scale=40).rvs(size=1000) # 1. Подгоняем нормальное и тяжелохвостое распределения Стьюдента к историческим данным norm_params = stats.norm.fit(latencies) t_params = stats.t.fit(latencies) # 2. Замораживаем подобранные модели fitted_norm = stats.norm(*norm_params) fitted_t = stats.t(*t_params) # 3. Вычисляем вероятность превышения порогового значения таймаута в 400 мс prob_outage_norm = fitted_norm.sf(400) prob_outage_t = fitted_t.sf(400) # 4. Рассчитайте время отклика 99-го процентиля для стресс-тестирования SLA. p99_norm = fitted_norm.ppf(0.99) p99_t = fitted_t.ppf(0.99) print(f»Нормальный SLA (99-й процентиль): {p99_norm:.2f} мс | P(Задержка > 400 мс): {prob_outage_norm:.4%}») print(f»SLA для тяжелых условий эксплуатации (99-й процентиль): {p99_t:.2f} мс | P(Задержка > 400 мс): {prob_outage_t:.4%}»)

Выход:

Нормальный SLA (99-й процентиль): 340,14 мс | P(Задержка > 400 мс): 0,0463% SLA при интенсивном t (99-й процентиль): 368,97 мс | P(Задержка > 400 мс): 0,6161%

Разница заметна. При наивном предположении о гауссовском распределении сбой с высокой задержкой, превышающей 400 мс, прогнозируется как редкое событие (происходящее в 0,0463% случаев). В действительности же модель с «тяжелыми хвостами» показывает, что вероятность сбоя составляет 0,6161% — более чем в 13 раз чаще. Подгонка распределений с «тяжелыми хвостами» предотвращает неожиданные сбои, которые кажутся «редкими».

# 5. Построение доверительных интервалов методом бутстреппинга для метрик сценария.

При проведении моделирования или анализе небольшого исторического набора данных вы рассчитываете суммарный показатель (например, 90-й процентиль размеров заказов или медианную себестоимость операций). Но насколько стабилен этот показатель? Что если ваша историческая выборка немного отличается?

Аналитический расчет доверительных интервалов для нестандартных показателей (таких как отношение или процентиль) является математически сложной задачей. Исторически сложилось так, что для ручной перевыборки данных специалисты использовали громоздкие вложенные циклы.

В SciPy 1.7 была представлена передовая функция scipy.stats.bootstrap. С помощью одного высокооптимизированного вызова функции можно вычислить надежные, непараметрические, скорректированные по смещению и ускоренные (BCa) доверительные интервалы для любой произвольной сводной статистики, без предположения о каком-либо базовом распределении.

Использование вложенных циклов для перевыборки, вычисления статистики и поиска квантилей бутстрап-распределения — медленный, многословный процесс, не позволяющий корректировать смещение или ускорение.

import numpy as np # Небольшая выборка из 50 сумм транзакций np.random.seed(42) transactions = np.random.lognormal(mean=4, sigma=0.8, size=50) # Цикл ручной бутстрап-проверки boot_medians = [] for _ in range(10000): sample = np.random.choice(transactions, size=len(transactions), replace=True) boot_medians.append(np.median(sample)) ci_low = np.percentile(boot_medians, 2.5) ci_high = np.percentile(boot_medians, 97.5) print(f»Доверительный интервал медианы ручной бутстрап-проверки: [{ci_low:.2f}, {ci_high:.2f}]»)

Выход:

Медианный доверительный интервал, полученный методом ручной бутстрап-методики: [36,47, 60,12]

В отличие от этого, передав наши данные и статистическую функцию непосредственно в stats.bootstrap, SciPy вычисляет скорректированный по смещению доверительный интервал за миллисекунды.

import numpy as np import scipy.stats as stats # Небольшая выборка из 50 сумм транзакций np.random.seed(42) transactions = np.random.lognormal(mean=4, sigma=0.8, size=50) # Определяем статистику, для которой хотим построить доверительный интервал (должна принимать данные в виде последовательности) def my_median(data_seq): return np.median(data_seq) # Выполняем stats.bootstrap с использованием метода BCa, передавая данные в виде кортежа, содержащего наш массив bootstrap_res = stats.bootstrap( (transactions,), statistic=my_median, n_resamples=9999, confidence_level=0.95, method='BCa' ) ci = bootstrap_res.confidence_interval print(f»Выборочный медианный размер транзакций: ${np.median(transactions):.2f}») print(f»95% доверительный интервал (BCa): [${ci.low:.2f}, ${ci.high:.2f}]») print(f»Стандартная ошибка медианы: ${bootstrap_res.standard_error:.4f}»)

Выход:

95% доверительный интервал (BCa): [$36,47, $60,12] Стандартная ошибка медианы: $5,8819

Обратите внимание, что метод BCa дал более узкий и высокоточный, математически обоснованный доверительный интервал по сравнению с простым ручным бутстрапом. BCa автоматически корректирует асимметрию и смещение в базовом наборе данных, обеспечивая принципиальную границу неопределенности для любого сценария или выборочной оценки.

# Завершение

Переход от простого мышления, основанного на точечных оценках, к зрелому дистрибутивному мышлению — один из самых эффективных шагов, которые может предпринять специалист по анализу данных.

Внедрив эти пять приемов из библиотеки scipy.stats в свой рабочий процесс моделирования, вы сможете проектировать высоконадежные, математически обоснованные системы сценариев:

• Замораживание распределений позволяет инкапсулировать ваши бизнес-предположения в чистые, взаимозаменяемые объекты сценариев.
• Метод Монте-Карло с использованием функции .rvs() обеспечивает беспрепятственное распространение многомерных неопределенностей в высокоскоростных векторизованных расширениях C.
• С помощью функции .ppf() выполняется параметрический анализ, позволяющий за микросекунды вычислить точные пороговые значения риска и процентили аналитически.
• Использование алгоритмов .fit() и .sf() для оптимизации с учетом «тяжелого хвоста распределения» защищает вашу производственную архитектуру от катастрофических событий типа «черный лебедь».
• Использование метода бутстрапа BCa с помощью stats.bootstrap позволяет получить надежные, не зависящие от распределения доверительные интервалы поверх любой метрики сценария.

В следующий раз, когда вам предстоит проводить стресс-тестирование систем или оценивать бизнес-риски в условиях неопределенности, забудьте о сложных внешних зависимостях. Возьмите стандартный набор инструментов Python для научных исследований, воспользуйтесь возможностями библиотеки scipy.stats и запустите моделирование!

Мэтью Мэйо ( @mattmayo13 ) имеет степень магистра компьютерных наук и диплом специалиста по анализу данных. Будучи главным редактором KDnuggets & Statology и внештатным редактором Machine Learning Mastery, Мэтью стремится сделать сложные концепции науки о данных доступными для всех. В сферу его профессиональных интересов входят обработка естественного языка, языковые модели, алгоритмы машинного обучения и изучение новых технологий искусственного интеллекта. Его движет стремление демократизировать знания в сообществе специалистов по науке о данных. Мэтью занимается программированием с 6 лет.

Источник: www.kdnuggets.com