python - 联合正态先验分布的后验
问题描述
我有一些关于高斯推理的基本问题。
我有以下数据:
(Log) dose, Number of animals, Number of deaths
-0.86, 5, 0
-0.30, 5, 1
-0.05, 5, 3
0.73, 5, 5
编辑:我假设剂量反应 logit(θ) = α + βx 的简单回归模型,其中 logit(θ) = log(θ / (1-θ))。θ 代表给定剂量 x 的死亡概率。
我想在 (α,β) 上创建一个联合正态先验分布,其中 α ∼ N(0,22),β ∼ N(10,102) 和 corr(α,β) = 0.5,然后计算 a 中的后验密度先验周围的点网格(α: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]])
这是正确的吗?
其次,如何计算网格中的后验密度?
解决方案
基于merv的好回答,回答我自己,我认为封闭的解决方案是:
p(yi|α,β,ni,xi)∝ [logit −1 (α+β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')
有一些差异。不确定分析优化函数是否正确。
希望这对其他人有帮助。
推荐阅读
- tensorflow - Tensorflow mIOU 和像素精度错误?
- azure - 如何使用 Powershell 获取与服务主体关联的证书的指纹?
- postgresql - 如何在 Pg Admin ||| 中运行 SQL 脚本而不丢失数据?
- r - 有没有办法在 R DT pacakge 中制作多级行分组表
- java - 通过无头安装 4J 自动化
- azure - 复制 blob 时 Azure SAS 到期
- inno-setup - 在存在重复项的 Inno Setup 中使用多个 DLL 文件
- android - 在路径上找不到类“com.jambod.user.model.ProcessModel”:仅在发布 apk 中出现 DexPathList 错误
- javascript - 为什么我无法在 ios9 门户网页上获取 localStorage?
- javascript - 响应式将文本输入内容右对齐