python - 在 PyMC3 中运行多元有序 logit

我正在尝试使用 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 使其运行？

最佳答案

在 之前我从未做过多变量序数回归，但似乎必须以两种方式解决建模问题:

预测器空间中的分区，在这种情况下，您需要切割线/曲线而不是点。

在已将预测变量空间投影到标量值的转换空间中进行分区，并且可以再次使用切割点。

如果你想使用 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

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

关于python - 在 PyMC3 中运行多元有序 logit，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/53198259/