Кзади от совместного нормального априорного распределения

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

p(yi|α,β,ni,xi)∝ [logit⁻¹(α+βxi)]^y * [1 − logit^{−1(α+β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')

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

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

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

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

Учитывая стандартные отклонения, 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)

Вычисление значений сетки или нет

Получение должным образом нормализованной апостериорной плотности требует маргинализации (интегрирования) по непрерывным переменным (например, θ), и это решаемо аналитически только в особых случаях, которые я не думаю, что ваш. Таким образом, вы можете либо вычислить интегралы и вычислить численные приближения, либо использовать какой-либо метод приближенного вывода, такой как 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 с помощью jitter+adapt_diag...

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

[beta] Выборка 2 цепочек: 100%|██████████| 220000/220000 [03:52‹00:00, 947,57 ничьих/с]

После настройки было 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()