在PyMC3中运行多元有序logit

Question

我试图用PyMC3建立一个贝叶斯多元有序logit模型。我根据this书中的例子得到了一个玩具多元逻辑模型。我还得到了一个基于this页面底部示例的有序逻辑回归模型。在

然而，我无法运行有序的多元逻辑回归。我认为问题可能在于切入点的指定方式，特别是形状参数，但我不确定如果存在多个独立变量，为什么会与只有一个独立变量不同，因为响应类别的数量没有改变。在

我的代码是：

MWE数据准备：

import pymc3 as pm
import numpy as np
import pandas as pd
from sklearn.datasets import load_iris

iris = load_iris(return_X_y=False)
iris = pd.DataFrame(data=np.c_[iris['data'], iris['target']],
                     columns=iris['feature_names'] + ['target'])

iris = iris.rename(index=str, columns={'sepal length (cm)': 'sepal_length', 'sepal width (cm)': 'sepal_width', 'target': 'species'})

以下是一个有效的多元（二进制）逻辑：

^{pr2}$

这是一个工作顺序逻辑图（带有一个自变量）：

x = iris['sepal_length'].values
y = pd.Categorical(iris['species']).codes

with pm.Model() as model:
    cutpoints = pm.Normal("cutpoints", mu=[-2,2], sd=10, shape=2,
                          transform=pm.distributions.transforms.ordered)

    y_ = pm.OrderedLogistic("y", cutpoints=cutpoints, eta=x, observed=y)
    tr = pm.sample(1000)

下面是我对多元有序logit的尝试，它打破了：

x = iris[['sepal_length', 'sepal_width']].values
y = pd.Categorical(iris['species']).codes

with pm.Model() as model:
    cutpoints = pm.Normal("cutpoints", mu=[-2,2], sd=10, shape=2,
                          transform=pm.distributions.transforms.ordered)

    y_ = pm.OrderedLogistic("y", cutpoints=cutpoints, eta=x, observed=y)
    tr = pm.sample(1000)

我得到的错误是：“ValueError：除了连接轴之外，所有输入数组的维数都必须完全匹配。”

这表明这是一个数据问题（x，y），但是数据看起来和multivariatorogit一样，这是有效的。在

如何修复有序多变量logit，使其运行？在

Answer 1

我以前从未做过多元序数回归，但似乎必须用两种方法来解决建模问题：

1. 在预测器空间中进行划分，在这种情况下，您需要剪切线/曲线，而不是点。在
2. 在已将预测值空间投影到标量值并可以再次使用切点的变换空间中的分区。在

如果您想使用pm.OrderedLogistic，那么似乎必须要使用后者，因为它似乎不支持开箱即用的多变量{}大小写。在

这是我的一个尝试，但我再次重申，我不确定这是一个标准的方法。在

import numpy as np
import pymc3 as pm
import pandas as pd
import theano.tensor as tt
from sklearn.datasets import load_iris

# Load data
iris = load_iris(return_X_y=False)
iris = pd.DataFrame(data=np.c_[iris['data'], iris['target']],
                     columns=iris['feature_names'] + ['target'])
iris = iris.rename(index=str, columns={
    'sepal length (cm)': 'sepal_length', 
    'sepal width (cm)': 'sepal_width', 
    'target': 'species'})

# Prep input data
Y = pd.Categorical(iris['species']).codes
X = iris[['sepal_length', 'sepal_width']].values

# augment X for simpler regression expression
X_aug = tt.concatenate((np.ones((X.shape[0], 1)), X), axis=1)

# Model with sampling
with pm.Model() as ordered_mvlogit:
    # regression coefficients
    beta = pm.Normal('beta', mu=0, sd=2, shape=X.shape[1] + 1)

    # transformed space (univariate real)
    eta = X_aug.dot(beta)

    # points for separating categories
    cutpoints = pm.Normal("cutpoints", mu=np.array([-1,1]), sd=1, shape=2,
                          transform=pm.distributions.transforms.ordered)

    y_ = pm.OrderedLogistic("y", cutpoints=cutpoints, eta=eta, observed=Y)

    trace_mvordlogit = pm.sample(5000)

这似乎收敛得很好，产生了不错的间隔

如果将beta和cutpoint的平均值返回到预测器空间，您将得到以下分区，这看起来是合理的。然而，萼片的长度和宽度并不是最好的分割。在

^{pr2}$

二阶条件

您可以很容易地扩展模型中的eta，以包括交互作用和高阶项，这样最终的分类器切割可以是曲线而不是简单的直线。例如，这里是二阶模型。在

from sklearn.preprocessing import scale

Y = pd.Categorical(iris['species']).codes

# scale X for better sampling
X = scale(iris[['sepal_length', 'sepal_width']].values)

# augment with intercept and second-order terms
X_aug = tt.concatenate((
    np.ones((X.shape[0], 1)), 
    X,
    (X[:,0]*X[:,0]).reshape((-1,1)),
    (X[:,1]*X[:,1]).reshape((-1,1)),
    (X[:,0]*X[:,1]).reshape((-1,1))), axis=1)

with pm.Model() as ordered_mvlogit_second:
    beta = pm.Normal('beta', mu=0, sd=2, shape=6)

    eta = X_aug.dot(beta)

    cutpoints = pm.Normal("cutpoints", mu=np.array([-1,1]), sd=1, shape=2,
                          transform=pm.distributions.transforms.ordered)

    y_ = pm.OrderedLogistic("y", cutpoints=cutpoints, eta=eta, observed=Y)

    trace_mvordlogit_second = pm.sample(tune=1000, draws=5000, chains=4, cores=4)

这个采样很好，所有的系数都有非零的hpd

如上所述，您可以生成分类区域的图