Гиббсовское сэмплирование на Python. Gibbs Sampling.. Gibbs Sampling. mcmc.. Gibbs Sampling. mcmc. python.. Gibbs Sampling. mcmc. python. цепи маркова.

Сегодня разбираем реализацию Gibbs Sampling на Python. Это один из методов Монте‑Карло по цепям Маркова (MCMC), который решает такую задачу:

«У нас есть сложное многомерное распределение, но мы не можем из него напрямую сэмплировать. Однако, если у нас есть условные распределения, то мы можем брать новые точки, обновляя поочередно каждую координату.»

Представим, что есть совместное распределение P(x, y). И неизвестно, как из него брать значения напрямую, но знаешь условные распределения:

P(x∣y),P(y∣x)

Как работает алгоритм Gibbs Sampling:

  • Выбираем начальные значения x(0),y(0)

  • Повторяем TTT итераций:

    • Генерируем x(t+1)∼P(x∣y(t))

    • Генерируем y(t+1)∼P(y∣x(t+1))

  • В итоге получаем выборку, которая сходится к истинному распределению P(x,y).

Алгоритм применяется в байесовском выводе, генеративных моделях, а также в некоторых финансовых моделях.

Как реализовать Gibbs Sampling в Python?

Представим, что есть совместное распределение:

P(x,y)∝exp(−(x 2 +y 2 +xy))

Где условные распределения:

P(x∣y)∼N(−0.5y,1)

P(y∣x)∼N(−0.5x,1)

Можно сэмплировать из этих распределений в Python.

Базовый Gibbs‑сэмплер:

import numpy as np
import matplotlib.pyplot as plt

# Определяем условные распределения
def sample_x_given_y(y):
    return np.random.normal(-0.5 * y, 1)

def sample_y_given_x(x):
    return np.random.normal(-0.5 * x, 1)

# Реализация Гиббсовского сэмплера
def gibbs_sampler(n_samples, init_x=0, init_y=0):
    x_samples = np.zeros(n_samples)
    y_samples = np.zeros(n_samples)

    x_samples[0] = init_x
    y_samples[0] = init_y

    for i in range(1, n_samples):
        x_samples[i] = sample_x_given_y(y_samples[i-1])
        y_samples[i] = sample_y_given_x(x_samples[i])

    return x_samples, y_samples

# Запускаем сэмплер
n_samples = 10000
x_samples, y_samples = gibbs_sampler(n_samples)

# Визуализируем результат
plt.scatter(x_samples, y_samples, alpha=0.2, s=5)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Гиббсовское сэмплирование")
plt.show()

Используем два условных распределения. Начинаем с какого‑то случайного (x,y)(x, y). По очереди обновляем xx и yy, используя условные распределения.

Выборка, которая имитирует истинное распределение P(x,y)

Выборка, которая имитирует истинное распределение P(x,y)

Пример прогноза цен на недвижимость

Вот мы работаем в компании, которая строит модель прогнозирования цен на квартиры. У нас есть исторические данные по цене квартир, площади, этажности, району и другим характеристикам.

Но проблема:

  • Данные шумные.

  • Линейная регрессия даёт ерунду (разброс слишком большой).

  • Хотим учесть неопределённость, чтобы строить доверительные интервалы.

Решение:

  • Байесовская линейная регрессия.

  • Оценка параметров с помощью Гиббсовского сэмплирования.

Математическая модель

Допустим, цена квартиры yi зависит от её площади x_i и района r_i:

yi​=α+βxi​+γri​+ϵi​,ϵi​∼N(0,σ2)

Где:

  • α— базовая цена.

  • β— цена за квадратный метр.

  • γ— коэффициент района (учитывает элитность).

  • σ2— дисперсия шума.

Мы не знаем параметры α,β,γ,σ, но можно оценить их с помощью Gibbs Sampling.

Генерация данных

Сгенерируем фейковые данные:

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# Генерируем 1000 квартир
n_samples = 1000

# Фактические параметры модели
alpha_true = 3_000_000  # Базовая цена квартиры
beta_true = 150_000     # Цена за м²
gamma_true = 500_000    # Разница в цене из-за района
sigma_true = 700_000    # Шум в цене

# Генерируем данные
x_area = np.random.uniform(30, 120, n_samples)  # Площадь от 30 до 120 м²
r_district = np.random.choice([0, 1], size=n_samples)  # Два района (0 = обычный, 1 = элитный)
epsilon = np.random.normal(0, sigma_true, n_samples)  # Шум

# Истинные цены квартир
y_price = alpha_true + beta_true * x_area + gamma_true * r_district + epsilon

# Визуализируем
plt.scatter(x_area, y_price, c=r_district, alpha=0.5)
plt.xlabel("Площадь (м²)")
plt.ylabel("Цена квартиры")
plt.title("Фейковые данные о ценах квартир")
plt.show()
Фейк.данные о ценах на квартиры

Фейк.данные о ценах на квартиры

Реализация Гиббсовского сэмплера

Теперь реализуем байесовскую регрессию, используя Гиббсовское сэмплирование.

def gibbs_sampler(n_iter, y, x, r):
    """Гиббсовский сэмплер для оценки параметров"""
    
    n = len(y)
    
    # Инициализация параметров
    alpha = 1_000_000
    beta = 100_000
    gamma = 200_000
    sigma = 1_000_000

    # Храним выборки параметров
    alpha_samples = np.zeros(n_iter)
    beta_samples = np.zeros(n_iter)
    gamma_samples = np.zeros(n_iter)
    sigma_samples = np.zeros(n_iter)

    for i in range(n_iter):
        # 1. Обновляем alpha
        alpha = np.random.normal(np.mean(y - beta*x - gamma*r), sigma / np.sqrt(n))

        # 2. Обновляем beta
        beta = np.random.normal(np.sum((y - alpha - gamma*r) * x) / np.sum(x**2), sigma / np.sqrt(np.sum(x**2)))

        # 3. Обновляем gamma
        gamma = np.random.normal(np.mean(y - alpha - beta*x), sigma / np.sqrt(n))

        # 4. Обновляем sigma
        residuals = y - (alpha + beta*x + gamma*r)
        sigma = np.sqrt(1 / np.random.gamma(n / 2, 2 / np.sum(residuals**2)))

        # Сохраняем результаты
        alpha_samples[i] = alpha
        beta_samples[i] = beta
        gamma_samples[i] = gamma
        sigma_samples[i] = sigma

    return alpha_samples, beta_samples, gamma_samples, sigma_samples

# Запускаем сэмплер
n_iter = 5000
alpha_samples, beta_samples, gamma_samples, sigma_samples = gibbs_sampler(n_iter, y_price, x_area, r_district)

Теперь посмотрим, как изменялись параметры:

plt.figure(figsize=(12, 6))
plt.subplot(2, 2, 1)
plt.plot(alpha_samples)
plt.title("Alpha (базовая цена)")

plt.subplot(2, 2, 2)
plt.plot(beta_samples)
plt.title("Beta (цена за м²)")

plt.subplot(2, 2, 3)
plt.plot(gamma_samples)
plt.title("Gamma (коэффициент района)")

plt.subplot(2, 2, 4)
plt.plot(sigma_samples)
plt.title("Sigma (шум)")

plt.tight_layout()
plt.show()

# Усредненные оценки параметров
alpha_est = np.mean(alpha_samples)
beta_est = np.mean(beta_samples)
gamma_est = np.mean(gamma_samples)
sigma_est = np.mean(sigma_samples)

print(f"Оцененные параметры:n alpha={alpha_est:.2f}, beta={beta_est:.2f}, gamma={gamma_est:.2f}, sigma={sigma_est:.2f}")
Динамика

Динамика

Графики показывает сходимость параметров в процессе Гиббсовского сэмплирования. В начале всех цепей наблюдается резкое изменение значений — это нормальное поведение из‑за плохого начального приближения. После ~1000 итераций параметры стабилизируются: Alpha (базовая цена) колеблется около 3 000 000, Beta (цена за м²) держится в районе 150 000, Sigma (шум) выходит на значение около 700 000, а Gamma (влияние района) остаётся самой нестабильной.

Нестабильность Gamma тут очевидно указывает на слабую информативность данных о районе или на недостаток наблюдений. Возможные решения: увеличить размер выборки, добавить новые признаки (например, количество комнат, год постройки) или пересмотреть саму модель. Визуально видно, что цепь сошлась, но стоит еще вдобавок увеличить burn‑in до 2000 итераций, чтобы убрать влияние первых нестабильных значений.

Так же можно подключить автокорреляцию цепей, построить гистограммы выборок параметров и сравнить их с истинными значениями. Если Gibbs‑сэмплирование работает медленно или нестабильно, стоит рассмотреть Hamiltonian Monte Carlo как альтернативу.

Хотите писать надежный код, который автоматически проверяет входные данные? Освойте Pydantic — мощную библиотеку для валидации и сериализации данных в Python. Познакомиться с Pydantic можно на открытом уроке 24 марта.

Список всех бесплатных уроков по ML, работе с данными и программированию смотрите в календаре мероприятий.

Автор: badcasedaily1

Источник

Версия для печати Версия для печати
Рейтинг@Mail.ru
Rambler's Top100