Я никогда не делал многомерную порядковую регрессию до , но кажется, что к проблеме моделирования нужно подходить двумя способами:
- Разделение в пространстве предиктора, в этом случае вам понадобятся линии / кривые пореза вместо точек.
- Разделение в преобразованном пространстве, в котором вы спроецировали пространство предиктора на скалярное значение и можете снова использовать точки отсечения.
Если вы хотите использовать 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
И, как указано выше, вы можете создать график классификации регионов
