python - 在 PyMC3 中运行多变量有序 logit

Question

我正在尝试使用 PyMC3 构建贝叶斯多变量有序 logit 模型。我已经根据 this 中的示例获得了一个玩具多元 logit 模型。书。我还得到了一个基于 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'})

这是一个有效的多变量(二进制)logit:

df = iris.loc[iris['species'].isin([0, 1])]
y = pd.Categorical(df['species']).codes
x = df[['sepal_length', 'sepal_width']].values

with pm.Model() as model_1:
      alpha = pm.Normal('alpha', mu=0, sd=10)
      beta = pm.Normal('beta', mu=0, sd=2, shape=x.shape[1])
      mu = alpha + pm.math.dot(x, beta)
      theta = 1 / (1 + pm.math.exp(-mu))
      y_ = pm.Bernoulli('yl', p=theta, observed=y)
      trace_1 = pm.sample(5000)

这是一个工作有序的 logit(带有一个自变量):

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)，但数据看起来与多元 logit 的数据相同，后者有效。

如何修复有序的多变量 logit 以便它运行？

Answer 1

我以前从未做过多元序数回归，但似乎必须以两种方式之一处理建模问题:

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

如果你想使用 pm.OrderedLogistic 似乎你必须使用后者，因为它似乎不支持多元 eta 案例盒子。

这是我的尝试，但同样，我不确定这是标准方法。

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 平均值返回到预测变量空间，您将得到以下分区，这看起来是合理的。然而，萼片的长度和宽度并不是最好的分区。

# Extract mean parameter values
b0, b1, b2 = trace_mvordlogit.get_values(varname='beta').mean(axis=0)
cut1, cut2 = trace_mvordlogit.get_values(varname='cutpoints').mean(axis=0)

# plotting parameters
x_min, y_min = X.min(axis=0)
x_max, y_max = X.max(axis=0)

buffer = 0.2
num_points = 37

# compute grid values
x = np.linspace(x_min - buffer, x_max + buffer, num_points)
y = np.linspace(y_min - buffer, y_max + buffer, num_points)

X_plt, Y_plt = np.meshgrid(x, y)
Z_plt = b0 + b1*X_plt + b2*Y_plt

# contour + scatter plots
plt.figure(figsize=(8,8))
plt.contourf(X_plt,Y_plt,Z_plt, levels=[-80, cut1, cut2, 50])
plt.scatter(iris.sepal_length, iris.sepal_width, c=iris.species)
plt.xlabel("Sepal Length")
plt.ylabel("Sepal Width")
plt.show()

二阶项

您可以轻松地扩展模型中的 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

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