Анализ сценариев с помощью распределений вероятностей SciPy
Монте-Карло, изображение предоставленоLuka NguyenотPixabay, бесплатно для коммерческого использования
В этом руководстве показано, как мы можем моделировать симуляцию Монте-Карло в Python. мы будем
- Используйте встроенные дистрибутивы SciPy, особенно. Нормальное, бета-распределение и распределение Вейбулла.
- Добавьте новый подкласс дистрибутива для дистрибутивов бета-PERT.
- Случайные числа получаются путем выборки латинского высокого куба.
- И построить три модели моделирования Монте-Карло.
0. Зависимость
Этот скрипт импортирует наши всепогодные пакеты pandas и numpy в дополнение к библиотеке SciPy _statistics_.
Концепция моделирования методом Монте-Карло
Изображение изPixabayизHans Braxmeier, бесплатно для коммерческого использования
Метод Монте-Карло (МК)— это метод моделирования, который строит распределения вероятностей для _выходных_ переменных модели, где некоторые _входные_ параметры являются случайными величинами. Метод MC иногда называют методом _многовероятностного моделирования_, потому что он включает в себя несколько случайных переменных, комбинированные эффекты которых нелегко описать одним уравнением в замкнутой форме.
Метод MC был изобретен Джоном фон Нейманом и Станиславом Уланом в конце 1940-х годов, когда они работали в лаборатории Лос-Аламоса. Когда они попытались рассчитать столкновения нейтронов с помощью детерминированного метода, они зашли в тупик. Идея Увама заключалась в том, чтобы использовать случайную выборку, основанную на раннем «суперкомпьютере» ENIAC, для получения приблизительных решений. Положения о правительственной тайне требовали кодового названия; фон Нейман и коллега Улана предложили название «Монте-Карло», где дядя Улана часто посещал казино перед Второй мировой войной: он «Просто поезжайте в Монте-Карло», чтобы сыграть в азартные игры.(Ułam)
Пример задачи моделирования MC: компания хочет смоделировать будущий успех нового продукта, поэтому она строит следующую имитационную модель.
- Ожидаемый объем продаж, v, является случайной величиной, которая соответствует бета-распределению PERT и оценивается по трем точкам.
- На цену продукта p будут влиять переговоры — некоторые клиенты будут просить скидки за количество или скидки за досрочную оплату; планировщики решают смоделировать неопределенность средней цены, которую получит фирма, как среднее значение 20 и стандартное отклонение 2. случайная переменная.
- На сырьевую стоимость продукта, m, повлияет предстоящее повышение цены поставщика, оцениваемое другим нормальным распределением.
- Модель должна обеспечивать несколько результатов: доход, прибыль, общие затраты.
Как правило, имитационная модель принимает входные параметры, определенные как случайные переменные, а затем возвращает распределение вероятностей для каждой выходной переменной, причем каждая выходная переменная является функцией своих входных параметров. Функция может быть простой суммой или произведением, экспоненциальным выражением входных данных или, как мы увидим в примере, функцией более высокого порядка, включающей вложенные случайные величины. Целевая переменная будет иметь функцию плотности, описывающую, как распределения входных переменных влияют друг на друга.
Имитационная модель выводит выходное распределение путем выборки входного распределения. 1000 или 100 000 итераций предоставили точки данных, описывающие поведение каждой входной переменной. Уравнение для выходной переменной преобразует эти входные данные в максимально возможное количество выборок выходной переменной. Получив этот образец, мы можем изучить его шаблон — его среднее значение, величину и склонность к выбросам — точно так же, как любой встроенный дистрибутив SciPy.
2. Сэмплирование латинского High Cube
Библиотека _random_ NumPy предоставляет методы для заполнения массивов одинаковыми случайными числами. Все объекты распределения SciPy имеют методы _rvs()_ для генерации случайных величин. Однако для моделирования методом Монте-Карло требуются большие массивы случайных чисел, для которых лучше всего подходят методы стратифицированной выборки, такие как выборка латинского гиперкуба (LHS).
LHS делит пространство выборки на непересекающиеся ячейки, каждая с равной вероятностью. LHS обеспечивает более равномерное извлечение случайных чисел из всего пространства выборки. Это предотвращает _скопление_ выбранных случайных чисел, которое может исказить частотную кривую. Выборка LHS приведет к более точным результатам моделирования. Это приведет к более низкой стандартной ошибке и уменьшит количество итераций.
Для демонстрации мы создаем небольшой массив из 10 случайных чисел с помощью сэмплера SciPy LatinHyperCube.
Равномерные (0;1) случайные числа после генерации могут масштабироваться для заполнения более широкого диапазона, например, от 0 до 100. Однако для наших целей будет достаточно стандартного случайного числа от 0 до 1.
В следующей главе мы будем использовать выборки LHS для генерации случайных величин с определенными распределениями: PERT, нормальные функции и функции Вейбулла.
3. Распределения вероятностей, которые моделируют неопределенность при моделировании
3.1 Распределение PERT для моделирования экспертных оценок
Моя предыдущая статья — «Моделирование экспертных оценок с помощью распределения Beta-PERT» — является объяснениемРаспределение PERTруководство(Анализ сценариев Python. Моделирование экспертных оценок с помощью распределения Beta-PERT | К науке о данных). )
Если эксперт в области, процессе или отрасли может предоставить так называемые _три балла_ для стохастического процессаоценка (три -point estimation), включая наихудший, возможный и наилучший результаты, функция PERT может справедливо соединить точки. Он преобразует эти точки в распределения вероятностей, которые мы можем назначить случайным переменным в модели моделирования.
Распределенный каталог SciPy 123 не включает функцию PERT. Поэтому мы создаем его как новый подкласс, наследуемый от родительского класса SciPy _rv_continuous_.
В нашем примере мы создаем экземпляр подкласса PERT с четырьмя параметрами для прогнозирования продаж нового продукта.
Мы выбрали это распределение PERT для имитации ожидаемых продаж с наиболее вероятным значением 12 000, заключенным между минимумом 8 000 и максимумом 18 000.
Его ожидаемая стоимость, в среднем, составит 12 333 проданных единицы. Оно имеет небольшую положительную асимметрию, все еще напоминающую нормальное распределение, и умеренно отрицательный избыточный эксцесс, равный -0,62, что делает его пластинчатым или громоздким распределением, склонным к большему количеству хвостов, чем нормальное распределение с небольшим количеством выбросов.
6 раундов по 60 на medium.com/Media/…
Строка 3 рисует выборку LHS из N = 10 000 однородных случайных величин. Затем строка 4 принимает массив в качестве входного параметра.процентный пунктФункция ppf()_ интерпретирует 10 000 стандартных средних чисел как вероятности, вычисляет 10 000 значений PERT и заносит их в массив randP.
На гистограмме показана форма смоделированного распределения PERT.
3.2 Нормальное распределение
Мы применяем тот же подход (за исключением необходимости создания нового подкласса распределения), беря 10 000 выборок LHS из нормального распределения, что, как мы надеемся, отражает неопределенность, присущую ценам реализации и затратам на единицу сырья. Средняя цена продажи составит 20 евро со стандартным отклонением 2 евро.
medium.com/Media/69 Фэйчанфен…
Второе нормальное распределение моделирует стоимость единицы сырья со средним значением 13 евро и стандартным отклонением 1,40 евро.
3.3 Моделирование 1 — Моделирование суммы и произведения случайных величин
Мы создаем три случайные величины, представляющие неопределенность трех факторов, которые будут определять вероятность того, что новый продукт принесет прибыль, или риск убытков.
- 1x PERT -- для продаж v.
- 2 нормальные переменные, представляющие продажную цену p и стоимость единицы сырья m.
Мы предполагаем, что другие затраты, не связанные с ценами поставщиков, могут быть отражены детерминированной переменной o.
Чтобы рассчитать целевую переменную имитационной модели — валовую прибыль продукта GP — мы связываем эти случайные переменные вместе следующим образом.
- v = количество, p = цена, m = стоимость единицы сырья, o = другая стоимость единицы продукции
- GP = v * (p - m - o)
В качестве вторичной выходной переменной мы можем моделировать GP одновременно с доходом R. В общем, как только мы определили их входные случайные переменные, мы можем настроить имитационные модели с любым количеством выходных переменных, которое мы хотим. Каждая выходная переменная может состоять из одной или нескольких входных случайных величин и будет описываться собственным массивом из 10000 переменных.
- R = v * p
Это устанавливает нашу простую имитационную модель.
Чтобы исследовать свойства массива, который возвращает распределение GP случайной величины валовой прибыли, мы написали функцию _dist_properties()_, которая будет считывать массив и возвращать его моменты и выбранные величины.
Для асимметрии и эксцесса мы используем функцию SciPy _skew().и эксцесс()_ нам нужно использовать метод преобразования _asscalar()_ numpy, чтобы преобразовать его результат из массива в плоское число. Собираем их значения в словарь _dict_moments_ и добавляем в него имя меры. Генератор списка в строке 14 выводит содержимое словаря построчно.
Затем мы вычисляем квантили ниже строки 17, собираем их в другой словарь _dict_quantiles_ и печатаем его через понимание списка в строке 27.
Мы объединяем эти два словаря, моменты и величины, чтобы сформировать всеобъемлющий словарь, _dict_metrics_.
medium.com/Media/147Чэнду…medium.com/Media/712ah4…
Средняя прибыль составит 49 250 евро, что близко к медиане.
И асимметрия, и эксцесс довольно умеренные, что означает, что по сравнению с нормальным распределением тенденция к выбросам незначительна.
Количества в таблице показывают, что если и цена продажи, и стоимость единицы сырья сходятся к наихудшему возможному результату, риск нового продукта составляет 5%. Прибыль превысит 10 731 евро с вероятностью 90% (размерность 10%).
В следующей главе мы пойдем еще дальше и построим имитационную модель более высокого порядка. В то время как первая модель собирала выходные переменные из сумм и произведений простых случайных величин, теперь мы разработаем имитационную модель, включающую _вложенные_ случайные переменные.
3.5 Моделирование 2. Вложенное распространение
Если вы читали мои предыдущие статьи(Intro to Probability Distributions and Distribution Fitting with Python's SciPy), вы помните, что мы применяли распределение Вейбулла для разрешения сбоеввремяилиСрок службы компонентапроблема. Срок службы компонента обычно можно использоватьWeibullМоделируется распределение, форма которого отражает _частоту отказов_, а параметр масштаба задает так называемый _характеристический ресурс_.
Первый космический корабль, который отправится к Марсу, будет оснащен электронными платами, срок службы которых может составлять 50 000 часов. Чтобы оценить резервирование, необходимое в конструкции космического корабля, мы смоделируем, какая часть всех цепей сгорит через сколько часов.
Предположим, что эквивалентно реальности, что два параметра распределения Вейбулла - форма и характерное время жизни - не являются детерминированными. Поэтому мы будем оценивать эти параметры распределения как наши собственные случайные величины. Форма Вейбулла и ее масштаб могут колебаться в пределах границ, установленных характером ее распределения.
- Команда инженеров подсчитала, что параметр формы должен иметь нормальное распределение со средним значением 1,5 и стандартным отклонением 0,1.
- Для характерного срока службы мы использовали оценку по трем точкам, предоставленную командой инженеров, максимальное, вероятное и минимальное значения при 45 000, 50 000 и 60 000 часов работы. Затем мы получаем распределение PERT, чтобы отразить диапазон неопределенности.
Каждая из двух вспомогательных функций отбирает по 10 000 случайных величин через латинские гиперкубы.
- _wei_shp()_ возвращает параметр формы нормального распределения.
- _wei_charlife()_ возвращает массив переменных PERT для времени жизни характеристики.
- Мы устанавливаем параметр положения Вейбулла (также известный как время ожидания) равным 0, что означает, что сбои продукта обнаруживаются, как только производство завершено и начинается проверка качества.
В строке 7 мы объединяем случайные величины в выходную переменную Вейбулла для расчета времени до отказа. Функция Вейбулла вызывает вспомогательные функции для получения параметров формы и масштаба. Целевая переменная _rand_CL_ будет содержать массив из 10 000 аналоговых выходов, представляющих срок службы компонента в часах.
Результаты с тремя вложенными распределениями: сумма PERT и нормальная случайная величина в качестве параметров распределения Вейбулла. Мы можем назвать этот результат гибким или случайным распределением Вейбулла, а не фиксированным распределением, параметры которого являются точечными значениями.
У нас есть два варианта просмотра свойств результатов.
- Проанализируйте его как массив, как мы делали с предыдущими результатами моделирования.
- Примените класс SciPy _rv_histogram, который делит выходной массив на гистограммы и превращает его в «истинное» распределение вероятностей SciPy, для которого мы можем вызывать функции распределения, такие как pdf и ppf.
Давайте взглянем на класс гистограммы.
На рисунке синим цветом показаны времена жизни бинов, которые мы смоделировали в массиве rand_CL. Оранжевым цветом показана функция плотности вероятности, полученная из массива с помощью класса _rv_histogram_.
Функция _histdist_properties()_ вычисляетсвойства распределений histdist_probability. Эта функция чем-то похожа на функцию _distribution_properties(), которую мы написали ранее.разные. Класс _rv_histogram_ не обеспечивает одни и те же свойства распределения во всех случаях. Например, минимальное и максимальное значения должны быть получены из метода ._support()_.
Мы вызываем _histdist_properties()функция для получения метрики в таблице. Фрейм данных содержит не только кумулятивную функцию распределения и величину, но и столбец случайных величин. функция выборки.rvs()_ позволяет нам выбрать несколько случайных чисел из нашего выходного распределения — при желании тысячи, но без стратифицированной выборки метода LHS.
Выше мы проанализировали непрерывную _rv_histogram_, полученную из массивараспределенный. Как видите, кривая распределения не совсем совпадает с данными массива.
Ниже мы вызываем функцию _dist_properties()_ для исходного массива из 10 000 выходных чисел и интерпретируем моменты и величины.
- Средний срок службы компонентов составит 45 647 часов.
- Естественно, продолжительность жизни будет смещена вправо с более чем 100 000 часов выбросов.
- Эксцесс уникален — распределение асимметрично со склонностью к долгоживущему компоненту в хвосте.
- Через 11 076 часов сгорит 10% компонентов.
Центр управления полетами сообщил нам, что 11 000 часов работы станут критическим порогом. После завершения этих часов использования цепь больше не будет критически важной. Эти цифры означают, что с точки зрения печатных плат космический корабль должен быть спроектирован с учетом как минимум 10-процентной избыточности: при добавлении запасных печатных плат они автоматически включаются взамен тех, которые начинают выходить из строя.
Итак, наш космический корабль почти готов к запуску.
фотоPeter HотPixabayСнято, бесплатно для коммерческого использования
3.6 Модель 3: Моделирование с фиксированными параметрами
В качестве последнего упражнения давайте рассмотрим, насколько модель Вейбулла чувствительна к неопределенностям в ее форме и сроке службы.
Мы запустим симуляцию Монте-Карло с распределением Вейбулла «фиксированной формы».
- Параметр формы будет зафиксирован на среднем значении 1,5.
- Время жизни функции: среднее значение распределения PERT. 50 500 часов.
Мы снова вызываем функцию _dist_properties()_ и вставляем ее метрики во фрейм данных, чтобы мы могли сравнить фиксированный Weibull рядом с ранее вложенным, поэтому «гибким» Weibull.
medium.com/Media/Storm становится…
Как и ожидалось, мы видим низкое стандартное отклонение для фиксированной переменной Вейбулла, которая была сгенерирована без неопределенности в параметрах формы и масштаба. Но разница небольшая. Ясно, что неопределенности в параметрах формы и масштаба гибкой модели Вейбулла довольно плотно сгруппированы вокруг среднего значения, которое мы используем в фиксированной модели Вейбулла. Меньшая асимметрия и эксцесс указывают на меньшую склонность к выбросам и асимметрии.
Моделирование с другими значениями параметров или разными распределениями может показать большие различия из-за повышенной неопределенности.
4.4 Заключение
На этом сегодняшний урок завершен.
Мы рассмотрели три примера симуляций Монте-Карло, созданных с помощью набора инструментов, предоставляемых библиотекой SciPy.
- Первая модель — имитация прибыли — показывает простые суммы и произведения случайных величин.
- Вторая симуляция — время до отказа Вейбулла — объясняет концепцию вложенных случайных величин, когда параметры распределения вероятностей сами являются случайными величинами.
- Третье моделирование делает смелое предположение: параметры распределения точно известны. Если мы игнорируем источники неопределенности в наших имитационных моделях, мы рискуем, однако, недооценить распределение возможных результатов, если только неопределенность не сведена к среднему значению параметров.
Блокноты Jupyter можно загрузить с GitHub:h3ik0th/МонтеКарлоСим. Моделирование методом Монте-Карло с помощью SciPy (github.com)
- Монте-Карло, изображение предоставленоLuka NguyenотPixabay, бесплатно для коммерческого использования
- Казино Монте-Карло, изображение предоставлено:Hans BraxmeierfromPixabay, бесплатно для коммерческого использования
- Фабрика, фото предоставлено:Peter HfromPixabay, бесплатно для коммерческого использования
- Все остальные изображения: предоставлены автором
Зависит отСимуляция Монте-Карло на основе PythonПервоначально опубликовано на MediumTowards Data Science, где люди продолжают беседу, выделяя историю и отвечая на нее.