Задний от сустава нормальное до распределения

Question

У меня есть несколько основных вопросов о гауссовском выводе.

У меня есть следующие данные:

(Log) dose, Number of animals, Number of deaths
-0.86, 5, 0
-0.30, 5, 1
-0.05, 5, 3
0.73, 5, 5

РЕДАКТИРОВАТЬ: Я предполагаю простую регрессионную модель для логарифмической зависимости доза-ответ (θ)= α + βx, где logit (θ) = log (θ / (1-θ)).θ обозначает вероятность смерти при данной дозе x.

Я хочу создать совместное нормальное предварительное распределение по (α, β) с α ∼ N (0,22), β ∼ N (10,102),и corr (α, β) = 0,5, а затем вычислите апостериорную плотность в сетке точек вокруг предыдущего (α: 0 ± 4, β: 10 ± 20).

Сначала я создал суставную нормальпредыдущее распределение следующее:

import numpy as np
from scipy import stats
x = np.array([-0.86, -0.30, -0.05, 0.73])
n = np.array([5, 5, 5, 5])
y = np.array([0, 1, 3, 5])
prior = stats.multivariate_normal([0, 10], [[0.5, 0], [0, 0.5]])

Это правильно?

Во-вторых, как рассчитать апостериорную плотность в сетке?

Answer 1

Основываясь на хорошем ответе Мерв, чтобы ответить самому себе, я думаю, что закрытое решение:

p (yi | α, β, ni, xi) ∝ [logit ⁻¹ (α + βxi)] ^y * [1 - logit ⁻¹ (α + βx) ^{n − y} ]

Таким образом, апостериор может быть рассчитан следующим образом:

import numpy as np
from scipy import optimize, stats
import matplotlib.pyplot as plt
x = np.array([-0.86, -0.30, -0.05, 0.73])
n = np.array([5, 5, 5, 5])
y = np.array([0, 1, 3, 5])

ngrid = 100
mu_1, mu_2, sd_1, sd_2 = 0, 10, 2**2, 10**2
A = np.linspace(-4, 4, ngrid)
B = np.linspace(-10, 30, ngrid)

mu = np.array([0, 10])
s = np.array([[22, 102]])
Rho = np.array([[1, 0.5], [0.5, 1]])
Sigma = Rho * np.outer(s, s)
prior = stats.multivariate_normal([mu_1, mu_2], Sigma)

def prop_likelihood(input_values):
    ilogit_abx = 1 / (np.exp(-(input_values[...,0][...,None]*np.ones(x.shape) + input_values[...,1][...,None] * x)) + 1)
    return np.prod(ilogit_abx**y * (1 - ilogit_abx)**(n - y), axis=ilogit_abx.ndim -1)

grid_a , grid_b = np.meshgrid(A,B)
grid = np.empty(grid_a.shape + (2,)) 
grid[:, :, 0] = grid_a
grid[:, :, 1] = grid_b

posterior_density = prior.pdf(grid)*prop_likelihood(grid)

Что тогда можно проиллюстрировать:

fig, ax = plt.subplots(figsize=(10, 5)
ax.imshow(
    posterior_density,
    origin='lower',
    aspect='auto',
    extent=(A[0], A[-1], B[0], B[-1])
)
ax.set_xlim([-4, 4])
ax.set_ylim([-10, 30])
ax.set_xlabel(r'$\alpha$')
ax.set_ylabel(r'$\beta$')
ax.set_title('Posterior heatmap')
ax.grid('off')

Аналитическое решение:

def opt(params):
    a, b  = params[0], params[1]
    z = np.exp(a + b * x) / (1 +  np.exp(a + b * x))
    e = - np.sum(y * np.log(z) + (n - y) * np.log(1 - z))
    return e

optim_res = optimize.minimize(opt, np.array([0.0, 0.0]))
mu_opt = optim_res['x']
sigma_opt = optim_res['hess_inv']
posterior_optimized = stats.multivariate_normal(mean=mu_opt, cov=sigma_opt)

Который затем можно построить

fig, ax = plt.subplots(figsize=(10, 5)
ax.imshow(
    posterior_optimized.pdf(grid),
    origin='lower',
    aspect='auto',
    extent=(A[0], A[-1], B[0], B[-1])
)
ax.set_xlim([-4, 4])
ax.set_ylim([-10, 30])
ax.set_xlabel(r'$\alpha$')
ax.set_ylabel(r'$\beta$')
ax.set_title('Posterior heatmap from analytical solution')
ax.grid('off')

Есть некоторые различия.Не уверен, что функция аналитической оптимизации верна.

Надеюсь, это поможет другим.

Answer 2

Параметризация Гаусса

Чтобы ответить на первый вопрос, вы неправильно параметрируете нормальное распределение.В частности, ваша ковариационная матрица не указана в соответствии с вашим описанием.

Учитывая стандартные отклонения s_1 = 22 и s_2 = 102 и желаемую корреляцию 0,5, правильная ковариационная матрица:

 ---                    ---
| s_1*s_1      0.5*s_1*s_2 |
|                          |
| 0.5*s_1*s_2      s_2*s_2 |
 ---                    ---

То есть отклонения по диагонали и ковариации вне диагонали.В Numpy / Scipy это будет

mu = np.array([0, 10])
s = np.array([[22, 102]])
Rho = np.array([[1, 0.5], [0.5, 1]])
Sigma = Rho * np.outer(s, s)

prior = stats.multivariate_normal(mean=mu, cov=Sigma)

---

Computing Grid Values ​​или Not

Для получения надлежащим образом нормализованной апостериорной плотности требуется маргинализация (интеграция) по непрерывным переменным (например, θ), и это аналитически разрешимо только в особых случаях, которые я не думаю, что у вас есть.Таким образом, вы можете либо вычислить интегралы и вычислить числовые приближения, либо использовать некоторый метод приближенного вывода, такой как MCMC или вариационный вывод.Для этого есть отличные инструменты, такие как PyMC3 и PyStan.

Получение апостериорных значений только для дискретных точек на сетке требует наложения условных значений на переменные модели.Однако большинство вероятностных инструментов программирования в наши дни настолько выразительны, что будет проще просто вывести полное апостериорное значение, а если у вас действительно есть какие-то особые значения сетки интереса, проверить их позже.

PyMC3 Пример

Вот полный задний вывод в PyMC3 с вашим сильным приоритетом:

import numpy as np
import pymc3 as pm
import theano
import theano.tensor as tt
import matplotlib.pyplot as plt
import arviz as az

# Data
X = np.array([-0.86, -0.30, -0.05, 0.73])
N = np.array([5, 5, 5, 5])
Y = np.array([0, 1, 3, 5])

# augment X for simpler regression expression
X_aug = tt.stack(np.ones_like(X), X).T

# Prior params
mu = np.array([0, 10])
sd = np.array([22, 102])
Rho = np.array([[1, 0.5],[0.5, 1]])
Sigma = np.outer(sd, sd) * Rho

with pm.Model() as binomial_regression:
    # regression coefficients (strong prior)
    beta = pm.MvNormal('beta', mu=mu, cov=Sigma, shape=2)

    # death probability
    theta_i = pm.Deterministic('theta_i', pm.invlogit(X_aug.dot(beta)))

    # outcomes
    y_i = pm.Binomial('y_i', n=N, p=theta_i, observed=Y)

    trace = pm.sample(10000, tune=100000, target_accept=0.8, random_seed=2018)

Это нормально, но выборка требует большого количества шагов настройки, чтобы уменьшить расхождения:

Автоматическое назначение сэмплера NUTS ...

Инициализация NUTS с использованием джиттера + adapt_diag ...

Многопроцессорная выборка (2 цепочки в 2 заданиях) NUTS:

[beta] Выборка из 2 цепей: 100% | ██████████ |220000/220000 [03:52 <00:00, 947,57 Draws / s] </p>

Было 1 расхождение после настройки.Увеличьте target_accept или выполните повторную параметризацию.

Число эффективных выборок для некоторых параметров меньше 25%.

Графики трассировки

Совместный участок

ax, _, _ = az.jointplot(trace, var_names=['beta'], kind='hexbin')
ax.set_xlabel("Intercept Coefficient ($\\beta_0$)")
ax.set_ylabel("Slope Coefficient ($\\beta_1$)")
plt.show()