Kodduu Python

Пример кода для проведения статистического анализа с использованием библиотеки PyStan (интерфейс Python для языка моделирования Stan), демонстрирующий выполнение байесовской линейной регрессии.


import pystan
import numpy as np
import matplotlib.pyplot as plt

# Генерация синтетических данных
np.random.seed(42)
size = 100
X = np.random.randn(size)
true_intercept = 1.0
true_slope = 2.5
true_sigma = 0.5
Y = true_intercept + true_slope * X + np.random.normal(0, true_sigma, size)

# Определяем модель на языке Stan
stan_model_code = """
data {
    int<lower=0> N;        // Количество наблюдений
    vector[N] x;           // Независимая переменная
    vector[N] y;           // Зависимая переменная
}
parameters {
    real alpha;            // Интерцепт
    real beta;             // Коэффициент наклона
    real<lower=0> sigma;   // Стандартное отклонение
}
model {
    y ~ normal(alpha + beta * x, sigma);  // Линейная модель с нормальным шумом
}
"""

# Подготовка данных для передачи в Stan
stan_data = {
    'N': size,
    'x': X,
    'y': Y
}

# Компиляция и запуск модели в Stan
stan_model = pystan.StanModel(model_code=stan_model_code)
fit = stan_model.sampling(data=stan_data, iter=2000, chains=4)

# Вывод результатов выборки
print(fit)

# Получение результатов выборки для параметров модели
fit_results = fit.extract()

# Визуализация апостериорных предсказаний
alpha_samples = fit_results['alpha']
beta_samples = fit_results['beta']

x_pred = np.linspace(-3, 3, 100)
y_pred = np.zeros((len(alpha_samples), len(x_pred)))

for i in range(len(alpha_samples)):
    y_pred[i, :] = alpha_samples[i] + beta_samples[i] * x_pred

# Построение графика с доверительными интервалами
plt.plot(x_pred, np.mean(y_pred, axis=0), label='Posterior mean', color='blue')
plt.fill_between(x_pred, np.percentile(y_pred, 2.5, axis=0), np.percentile(y_pred, 97.5, axis=0), color='blue', alpha=0.3)
plt.scatter(X, Y, c='r', label='Observed data')
plt.title('Байесовская линейная регрессия с PyStan')
plt.legend()
plt.show()

### Описание анализа

#### Шаги:
1. Генерация данных: Мы создаем синтетические данные с линейной зависимостью между независимой переменной \(X\) и зависимой переменной \(Y\) с добавлением случайного шума.
2. Модель на языке Stan: Мы описываем байесовскую линейную регрессию с нормальным распределением шума в модели. В этой модели есть три параметра: интерцепт (alpha), коэффициент наклона (beta) и стандартное отклонение ошибки (sigma).
3. Подготовка данных для Stan: Мы передаем данные, такие как количество наблюдений и сами данные, в формате, понятном Stan.
4. Запуск выборки: Мы используем алгоритм HMC (Гамильтонова Монте-Карло) для выполнения байесовского вывода с использованием выборки параметров модели. Параметры выборки будут использованы для построения апостериорных распределений.
5. Визуализация результатов: Мы визуализируем предсказания модели, построенные на основе апостериорных распределений, с доверительными интервалами. Это позволяет оценить неопределенность предсказаний.

### Особенности PyStan:
1. Высокая производительность: PyStan является оберткой над Stan, который использует высокоэффективные алгоритмы Гамильтоновой Монте-Карло для проведения выборки из апостериорных распределений. Это делает его особенно полезным для сложных и вычислительно тяжелых моделей.
2. Явное описание модели: Stan использует свой язык для определения моделей, что делает описания вероятностных моделей более строгими и легко читаемыми. Это позволяет четко задавать как априорные распределения, так и зависимость между переменными.
3. Поддержка сложных иерархических моделей: PyStan может использоваться для моделирования сложных иерархических моделей, что делает его подходящим для задач с большим количеством параметров и зависимостей.
4. Мощные методы байесовского вывода: Используя такие алгоритмы, как NUTS и HMC, PyStan обеспечивает эффективный байесовский анализ, особенно в пространствах параметров высокой размерности.

Подпишись 👉🏻 @KodduuPython 🤖

🔥2

162 views07:10