PyMC3: байесовское моделирование и прогнозирование в Python

Привет, Хабр!

Сегодня мы рассмотрим то, как реализовать байесовское моделирование и прогнозирование с использованием замечательной библиотеки PyMC3.

Байесовские методы — подход к статистическому моделированию, который включает в себя оценку вероятностных распределений параметров модели на основе данных.

Установим

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

pip install pymc3

conda install -c conda-forge pymc3

Для полноценной работы PyMC3 требуется несколько зависимостей:

1. Theano: основа для вычислений

   pip install Theano

2. ArviZ: для анализа и визуализации байесовских моделей

   pip install arviz

3. Дополнительные пакеты: NumPy, SciPy, которые являются стандартными библиотеками для работы с данными в Python.

   pip install numpy scipy

После установки всего можем приступать.

Создание байесовских моделей с PyMC3

Простая линейная регрессия

Линейная регрессия — это один из самых простых и распространенных методов статистического моделирования. В байесовском подходе мы задаем априорные распределения для параметров модели и используем данные для обновления наших знаний о параметрах.

Пример кода с PyMC3:

import pymc3 as pm import numpy as np import matplotlib.pyplot as plt  # генерация данных np.random.seed(42) X = np.random.randn(100) Y = 3 * X + np.random.randn(100)  # построение модели with pm.Model() as model:     # априорные распределения     alpha = pm.Normal('alpha', mu=0, sigma=10)     beta = pm.Normal('beta', mu=0, sigma=10)     sigma = pm.HalfNormal('sigma', sigma=1)      # линейная зависимость     mu = alpha + beta * X      # наблюдения     Y_obs = pm.Normal('Y_obs', mu=mu, sigma=sigma, observed=Y)      # монтекарловское семплирование     trace = pm.sample(2000, return_inferencedata=False)  # визуализация  pm.traceplot(trace) plt.show()

Задаем априорные распределения для параметров alpha, beta и sigma. Затем определяем линейную зависимость mu, и связываем наблюдаемые данные Y с помощью наблюдаемого распределения Y_obs.

Моделирование временных рядов

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

Пример с PyMC3:

import pymc3 as pm import numpy as np import matplotlib.pyplot as plt  # генерация данных временного ряда np.random.seed(42) X = np.linspace(0, 10, 100) Y = np.sin(X) + np.random.randn(100) * 0.1  # построение модели with pm.Model() as model:     # определение ядра гауссовского процесса     cov = pm.gp.cov.ExpQuad(1, ls=1)          # определение GP     gp = pm.gp.Marginal(cov_func=cov)          # наблюдения     sigma = pm.HalfNormal("sigma", sigma=0.5)     y_ = gp.marginal_likelihood("y", X=X[:, None], y=Y, noise=sigma)          # монтекарловское семплирование     trace = pm.sample(1000, return_inferencedata=False)  # визуализация with model:     y_pred = gp.conditional("y_pred", X[:, None])     samples = pm.sample_posterior_predictive(trace, vars=[y_pred], samples=200)  plt.plot(X, Y, 'ok', ms=3, alpha=0.5, label="Observed data") plt.plot(X, np.mean(samples['y_pred'], axis=0), label="Mean prediction") plt.fill_between(X, np.percentile(samples['y_pred'], 2.5, axis=0), np.percentile(samples['y_pred'], 97.5, axis=0), alpha=0.3, label="95% credible interval") plt.legend() plt.show()

Создали гауссовский процесс с использованием экспоненциально-квадратичного ядра и используем его для моделирования и прогнозирования временного ряда.

Иерархические модели

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

Пример:

import pymc3 as pm import numpy as np  # генерация данных np.random.seed(42) n_groups = 10 n_samples = 100 groups = np.random.randint(0, n_groups, n_samples) X = np.random.randn(n_samples) Y = 3 * X + groups + np.random.randn(n_samples)  # построение модели with pm.Model() as model:     # априорные распределения на уровне групп     group_mu = pm.Normal('group_mu', mu=0, sigma=10, shape=n_groups)     group_sigma = pm.HalfNormal('group_sigma', sigma=1)          # априорные распределения на глобальном уровне     alpha = pm.Normal('alpha', mu=0, sigma=10)     beta = pm.Normal('beta', mu=0, sigma=10)     sigma = pm.HalfNormal('sigma', sigma=1)          # линейная зависимость     mu = alpha + beta * X + group_mu[groups]          # наблюдения     Y_нobs = pm.Normal('Y_obs', mu=mu, sigma=sigma, observed=Y)          # монтекарловское семплирование     trace = pm.sample(2000, return_inferencedata=False)

group_mu моделирует случайные эффекты на уровне групп, а alpha и beta — глобальные эффекты.

Проверка и оценка моделей

В PyMC3 есть инструменты для диагностики и оценки моделей: forest plots и trace plots.

Пример trace plots:

import pymc3 as pm import matplotlib.pyplot as plt  # визуализация трассы pm.traceplot(trace) plt.show()

Пример forest plots:

import arviz as az  # визуализация forest plot az.plot_forest(trace, var_names=["alpha", "beta"], combined=True) plt.show()

Подробнее с PyMC3 можно ознакомиться здесь.

Больше про инструменты аналитики эксперты OTUS рассказывают в рамках практических онлайн-курсов. С полным каталогом курсов можно ознакомиться по ссылке.