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

Сегодня мы рассмотрим то, как реализовать байесовское моделирование и прогнозирование с использованием замечательной библиотеки 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 рассказывают в рамках практических онлайн-курсов. С полным каталогом курсов можно ознакомиться по ссылке.