Привет, Хабр!
Сегодня мы рассмотрим то, как реализовать байесовское моделирование и прогнозирование с использованием замечательной библиотеки 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 рассказывают в рамках практических онлайн-курсов. С полным каталогом курсов можно ознакомиться по ссылке.
ссылка на оригинал статьи https://habr.com/ru/articles/829874/
Добавить комментарий