0

我尝试使用 PyMC 解决逻辑回归模型。但是，诊断图显示了非常高的自相关性，并且在从后验分布重复采样后，我有时会得到非常不同的结果，所以我可能没有正确使用 PyMC。

模型如下：

Y_i = Bernoulli(p_i)
logit(p_i) = b0 + b1*x1 + b2*x2

实现是这样的：

import pymc as pm
import numpy as np

b0 =  pm.Normal('b0',  0., 1e-6, value=0.)
b1 =  pm.Normal('b1',  0., 1e-6, value=0.)
b2 =  pm.Normal('b2',  0., 1e-6, value=0.)

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32])
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15])
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0])

@pm.deterministic
def p(b0=b0, b1=b1, b2=b2, x1=x1, x2=x2):
    return np.exp(b0 + b1*x1 + b2*x2)/(1.0 + np.exp(b0 + b1*x1 + b2*x2))

obs = pm.Bernoulli('obs', p=p, value=value, observed=True)

m = pm.MCMC([obs, b0, b1, b2])

当我采样m.sample(500000, 200000, 50)并绘制结果后验分布时，我得到了：

在此处输入图像描述 在此处输入图像描述 在此处输入图像描述

在第二次尝试获得更好的结果时，我使用了@pm.observed：

import pymc as pm
import numpy as np

b0 =  pm.Normal('b0',  0., 1e-6, value=0.)
b1 =  pm.Normal('b1',  0., 1e-6, value=0.)
b2 =  pm.Normal('b2',  0., 1e-6, value=0.)

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32])
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15])
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0])

p = b0 + b1*x1 + b2*x2

@pm.observed
def obs(p=p, value=value):
    return pm.bernoulli_like(value, pm.invlogit(p))

m = pm.MCMC([obs, p, b0, b1, b2])

但它也会产生高自相关。

我增加了样本量但没有取得多大成功。我错过了什么？

4

1 回答 1

2

你肯定有一些收敛问题。部分问题是您的样本量非常小（n = 9），但正在尝试拟合 3 个参数。通过块更新参数，我获得了更好的里程数，但截距仍然很差：

beta =  pm.Normal('beta',  0., 0.001, value=[0.]*3)

x1 = np.array([31, 31, 36, 30, 32, 33, 31, 33, 32])
x2 = np.array([31, 16, 35, 42, 19, 37, 29, 23, 15])
value = np.array([1, 1, 1, 1, 0, 1, 1, 1, 0])

p = pm.Lambda('p', lambda b=beta: pm.invlogit(b[0] + b[1]*x1 + b[2]*x2))

obs = pm.Bernoulli('obs', p=p, value=value, observed=True)

（顺便说一句，当您对参数进行 logit 转换时，您的 beta 不需要超扩散先验——小于 0.001 是毫无结果的）。通过在 beta 上使用自适应 Metropolis 采样器，事情得到了进一步的改进：

M.use_step_method(AdaptiveMetropolis, M.beta)

这会导致收敛，但如您所见，没有信息可以告知先验：

概括

于 2014-09-01T20:09:43.033 回答