PyMC3: байесовское моделирование и прогнозирование в Python
![3fc8f0aaa3824555cf497d387e9ebe48.png](https://habrastorage.org/getpro/habr/upload_files/3fc/8f0/aaa/3fc8f0aaa3824555cf497d387e9ebe48.png)
Привет, Хабр!
Сегодня мы рассмотрим то, как реализовать байесовское моделирование и прогнозирование с использованием замечательной библиотеки PyMC3.
Байесовские методы — подход к статистическому моделированию, который включает в себя оценку вероятностных распределений параметров модели на основе данных.
Установим
Первый шаг к использованию PyMC3 — это, конечно же, установка библиотеки. PyMC3 можно установить двумя способами: через pip
и conda
:
pip install pymc3
conda install -c conda-forge pymc3
Для полноценной работы PyMC3 требуется несколько зависимостей:
Theano: основа для вычислений
pip install Theano
ArviZ: для анализа и визуализации байесовских моделей
pip install arviz
Дополнительные пакеты: 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 рассказывают в рамках практических онлайн-курсов. С полным каталогом курсов можно ознакомиться по ссылке.