有序多分类logistic模型python案例 有序多分类回归

目录：

• 第3章 线性分类

• 3.2 基于Softmax回归的多分类任务

• 3.2.1 数据集构建
• 3.2.2 模型构建

• 3.2.2.1 Softmax函数
• 3.2.2.2 Softmax回归算子

• 3.2.3 损失函数
• 3.2.4 模型优化

• 3.2.4.1 梯度计算
• 3.2.4.2 参数更新

• 3.2.5 模型训练
• 3.2.6 模型评价

第3章 线性分类

注： 这篇内容接上一篇，该篇代码有的调用了上篇的函数。

3.2 基于Softmax回归的多分类任务

Logistic回归可以有效地解决二分类问题，但在分类任务中，还有一类多分类问题，即类别数C大于2 的分类问题。Softmax回归就是Logistic回归在多分类问题上的推广。

使用Softmax回归模型对一个简单的数据集进行多分类实验。

3.2.1 数据集构建

我们首先构建一个简单的多分类任务，并构建训练集、验证集和测试集。
本任务的数据来自3个不同的簇，每个簇对一个类别。我们采集1000条样本，每个样本包含2个特征。

数据集的构建函数make_multi的代码实现如下：

import numpy as np

# 数据集的构建函数make_multi的代码实现如下：
def make_multiclass_classification(n_samples=100, n_features=2, n_classes=3, shuffle=True, noise=0.1):
    '''
     生成带噪音的多类别数据
    输入：
        - n_samples：数据量大小，数据类型为int
        - n_features：特征数量，数据类型为int
        - shuffle：是否打乱数据，数据类型为bool
        - noise：以多大的程度增加噪声，数据类型为None或float，noise为None时表示不增加噪声
    输出：
        - X：特征数据，shape=[n_samples,2]
        - y：标签数据, shape=[n_samples,1]
    '''
    # 计算每个类别的样本数量
    n_samples_per_class = [int(n_samples / n_classes) for k in range(n_classes)]
    for i in range(n_samples - sum(n_samples_per_class)):
        n_samples_per_class[i % n_classes] += 1
    # 将特征和标签初始化为0
    X = torch.zeros([n_samples, n_features])
    y = torch.zeros([n_samples], dtype = torch.int32)
    # 随机生成3个簇中心作为类别中心
    centroids = torch.randperm(2 ** n_features)[:n_classes]
    centroids_bin = np.unpackbits(centroids.numpy().astype('uint8')).reshape((-1, 8))[:, -n_features:]
    centroids = torch.tensor(centroids_bin, dtype=torch.float32)
    # 控制簇中心的分离程度
    centroids = 1.5 * centroids - 1
    # 随机生成特征值
    X[:, :n_features] = torch.randn(size=[n_samples, n_features])

    stop = 0
    for k, centroid in enumerate(centroids):
        start, stop = stop, stop + n_samples_per_class[k]
        # 指定标签值
        y[start:stop] = k % n_classes
        X_k = X[start:stop, :n_features]
        # 控制每个类别特征值的分散程度
        A = 2 * torch.rand(size=[n_features, n_features]) - 1
        X_k[...] = torch.matmul(X_k, A)
        X_k += centroid
        X[start:stop, :n_features] = X_k

    # 如果noise不为None，则给特征加入噪声
    if noise > 0.0:
        # 生成noise掩膜，用来指定给那些样本加入噪声,随机数据，< noise加噪音 
        noise_mask = torch.rand([n_samples]) < noise 
        for i in range(len(noise_mask)):
            if noise_mask[i]:
                # 给加噪声的样本随机赋标签值
                y[i] = torch.randint(n_classes, size=[1]).int()
    # 如果shuffle为True，将所有数据打乱
    if shuffle:
        idx = torch.randperm(X.shape[0])
        X = X[idx]
        y = y[idx]

    return X, y

随机采集1000个样本，并进行可视化。

# 固定随机种子，保持每次运行结果一致
torch.manual_seed(102)
# 采样1000个样本
n_samples = 1000
X, y = make_multiclass_classification(n_samples=n_samples, n_features=2, n_classes=3, noise=0.2)

# 可视化生产的数据集，不同颜色代表不同类别
plt.figure(figsize=(5,5))
plt.scatter(x=X[:, 0].tolist(), y=X[:, 1].tolist(), marker='*', c=y.tolist())
plt.savefig('linear-dataset-vis2.pdf')
plt.show()

运行结果：

将实验数据拆分成训练集、验证集和测试集。其中训练集640条、验证集160条、测试集200条。

# 将实验数据拆分成训练集、验证集和测试集。其中训练集640条、验证集160条、测试集200条。
num_train = 640
num_dev = 160
num_test = 200

X_train, y_train = X[:num_train], y[:num_train]
X_dev, y_dev = X[num_train:num_train + num_dev], y[num_train:num_train + num_dev]
X_test, y_test = X[num_train + num_dev:], y[num_train + num_dev:]

# 打印X_train和y_train的维度
print("X_train shape: ", X_train.shape, "y_train shape: ", y_train.shape)

运行结果：

X_train shape:  torch.Size([640, 2]) y_train shape:  torch.Size([640])

这样，我们就完成了Multi1000数据集的构建。

# 打印前5个数据的标签
print(y_train[:5])

运行结果：

tensor([0, 2, 0, 0, 1], dtype=torch.int32)

3.2.2 模型构建

在Softmax回归中，对类别进行预测的方式是预测输入属于每个类别的条件概率。与Logistic 回归不同的是，Softmax回归的输出值个数等于类别数 $C$，而每个类别的概率值则通过Softmax函数进行求解。

3.2.2.1 Softmax函数

Softmax函数可以将多个标量映射为一个概率分布。对于一个$K$维向量，$x=[x_1,⋯,x_K]$，Softmax的计算公式为：
$softmax(x_k)={exp(xk)\over ∑^K_{i=1}exp(xi)}。（3.11）$在Softmax函数的计算过程中，要注意上溢出和下溢出的问题。假设Softmax 函数中所有的 $x_k$ 都是相同大小的数值$a$，理论上，所有的输出都应该为$1\over k$。但需要考虑如下两种特殊情况：

• $a$为一个非常大的负数，此时exp(a) 会发生下溢出现象。计算机在进行数值计算时，当数值过小，会被四舍五入为0。此时，Softmax函数的分母会变为0，导致计算出现问题；
• $a$为一个非常大的正数，此时会导致exp(a)发生上溢出现象，导致计算出现问题。

为了解决上溢出和下溢出的问题，在计算Softmax函数时，可以使用$x_k−max(x)$代替$x_k$。 此时，通过减去最大值，$x_k$

Softmax函数的代码实现如下：

def softmax(X):
    """
    输入：
        - X：shape=[N, C]，N为向量数量，C为向量维度
    """
    x_max = torch.max(X, dim=1, keepdim=True).values # N,1
    x_exp = torch.exp(X - x_max)
    partition = torch.sum(x_exp, dim=1, keepdim=True) # N,1
    return x_exp / partition

# 观察softmax的计算方式
X = torch.tensor([[0.1, 0.2, 0.3, 0.4],[1,2,3,4]])
predict = softmax(X)
print(predict)

运行结果：

tensor([[0.2138, 0.2363, 0.2612, 0.2887],
        [0.0321, 0.0871, 0.2369, 0.6439]])

3.2.2.2 Softmax回归算子

在Softmax回归中，类别标签 $y∈\{1,2,…,C\}$ 。给定一个样本 $x$ ，使用Softmax回归预测的属于类别 $c$ 的条件概率为
$p(y=c|x)=softmax(w^T_cx+b_c)，（3.12）$其中$w_c$是第 $c$ 类的权重向量，$b_c$是第 $c$

Softmax回归模型其实就是线性函数与Softmax函数的组合。

将N个样本归为一组进行成批地预测。

$\hat Y=softmax(XW+b),(3.13)$其中$X∈R^{N×D}$为 $N$ 个样本的特征矩阵，$W=[w_1,……,w_C]$为 $C$ 个类的权重向量组成的矩阵，$\hat Y∈R^C$为所有类别的预测条件概率组成的矩阵。

我们根据公式（3.13）实现Softmax回归算子，代码实现如下：

class model_SR(op.Op):
    def __init__(self, input_dim, output_dim):
        super(model_SR, self).__init__()
        self.params = {}
        # 将线性层的权重参数全部初始化为0
        self.params['W'] = torch.zeros(size=[input_dim, output_dim])
        # self.params['W'] = paddle.normal(mean=0, std=0.01, shape=[input_dim, output_dim])
        # 将线性层的偏置参数初始化为0
        self.params['b'] = torch.zeros(size=[output_dim])
        self.outputs = None
    
    def __call__(self, inputs):
        return self.forward(inputs)

    def forward(self, inputs):
        """
        输入：
            - inputs: shape=[N,D], N是样本数量，D是特征维度
        输出：
            - outputs：预测值，shape=[N,C]，C是类别数
        """
        # 线性计算
        score = torch.matmul(inputs, self.params['W']) + self.params['b']
        # Softmax函数
        self.outputs = softmax(score)
        return self.outputs

# 随机生成1条长度为4的数据
inputs = torch.randn(size=[1, 4])
print('Input is:', inputs)
# 实例化模型，这里令输入长度为4，输出类别数为3
model = model_SR(input_dim=4, output_dim=3)
outputs = model(inputs)
print('Output is:', outputs)

运行结果：

Input is: tensor([[ 1.2374, -0.4203,  1.5545,  0.2168]])
Output is: tensor([[0.3333, 0.3333, 0.3333]])

从输出结果可以看出，采用全0初始化后，属于每个类别的条件概率均为$1\over C$。这是因为，不论输入值的大小为多少，线性函数 $f(x;W,b)$

3.2.3 损失函数

Softmax回归同样使用交叉熵损失作为损失函数，并使用梯度下降法对参数进行优化。通常使用 $C$ 维的one-hot类型向量$y∈\{0,1\}^C$ 来表示多分类任务中的类别标签。对于类别 $c$ ，其向量表示为：
$y=[I(1=c),I(2=c),…,I(C=c)]^T，（3.14）$其中$I(⋅)$是指示函数，即括号内的输入为“真”，$I(⋅)=1$；否则，$I(⋅)=0$。

给定有N个训练样本的训练集${(x^{(n)},y^{(n)})}^N_{n=1}$，令$\hat y^{(n)}=softmax(W^Tx^{(n)}+b)$为样本$x^{(n)}$在每个类别的后验概率。

多分类问题的交叉熵损失函数定义为：
$R(W,b)=−{1\over N}∑_{n=1}^N(y^{(n)})^Tlog\hat y^{(n)}=−{1\over N}∑_{n=1}^N∑_{c=1}^Cy^{(n)}clog \hat y^{(n)}_c.（3.15）$观察上式，$y^{(n)}_c$在c为真实类别时为1，其余都为0。也就是说，交叉熵损失只关心正确类别的预测概率，因此，上式又可以优化为：
$R(W,b)=−{1\over N}∑_{n=1}^Nlog[\hat y^{(n)}]_{y^{(n)}},（3.16）$其中$y^{(n)}$是第n个样本的标签。

因此，多类交叉熵损失函数的代码实现如下：

# 因此，多类交叉熵损失函数的代码实现如下:
class MultiCrossEntropyLoss(op.Op):
    def __init__(self):
        self.predicts = None
        self.labels = None
        self.num = None
    
    def __call__(self, predicts, labels):
        return self.forward(predicts, labels)

    def forward(self, predicts, labels):
        """
        输入：
            - predicts：预测值，shape=[N, 1]，N为样本数量
            - labels：真实标签，shape=[N, 1]
        输出：
            - 损失值：shape=[1]
        """
        self.predicts = predicts
        self.labels = labels
        self.num = self.predicts.shape[0]
        loss = 0
        for i in range(0, self.num):
            index = self.labels[i]
            loss -= torch.log(self.predicts[i][index])
        return loss / self.num

# 测试一下
# 假设真实标签为第1类
labels = torch.tensor([0])
# 计算风险函数
mce_loss = MultiCrossEntropyLoss()
print(mce_loss(outputs, labels))

运行结果：

tensor(1.0986)

3.2.4 模型优化

使用梯度下降法进行参数学习。

3.2.4.1 梯度计算

计算风险函数 $R(W,b)$ 关于参数 $W$ 和 $b$ 的偏导数。在Softmax回归中，计算方法为：
${∂R(W,b)\over ∂W}=−{1\over N}∑_{n=1}^Nx^{(n)}(y^{(n)}−\hat y^{(n)})^T=−{1\over N}X^T(y−\hat y),（3.17）$${∂R(W,b) \over ∂b}=−{1\over N}∑_{n=1}^N(y^{(n)}−\hat y^{(n)})^T=−{1\over N}1(y−\hat y).（3.18）$其中$X∈R^{N×D}$为N个样本组成的矩阵，$y∈R^N$ 为$N$个样本标签组成的向量，$\hat y∈R^N$为 $N$ 个样本的预测标签组成的向量，1为 $N$

将上述计算方法定义在模型的backward函数中，代码实现如下：

class model_SR(op.Op):
    def __init__(self, input_dim, output_dim):
        super(model_SR, self).__init__()
        self.params = {}
        # 将线性层的权重参数全部初始化为0
        self.params['W'] = torch.zeros(size=[input_dim, output_dim])
        # self.params['W'] = paddle.normal(mean=0, std=0.01, shape=[input_dim, output_dim])
        # 将线性层的偏置参数初始化为0
        self.params['b'] = torch.zeros(size=[output_dim])
        # 存放参数的梯度
        self.grads = {}
        self.X = None
        self.outputs = None
        self.output_dim = output_dim

    def __call__(self, inputs):
        return self.forward(inputs)

    def forward(self, inputs):
        self.X = inputs
        # 线性算子
        score = torch.matmul(self.X, self.params['W']) + self.params['b']
        # Softmax函数
        self.outputs = softmax(score)
        return self.outputs

    def backward(self, labels):
        """
        输入：
            - labels：真实标签，shape=[N, 1]，其中N为样本数量
        """
        # 计算偏导数
        N =labels.shape[0]
        labels = torch.nn.functional.one_hot(labels.to(torch.int64), self.output_dim)
        self.grads['W'] = -1 / N * torch.matmul(self.X.t(), (labels-self.outputs))
        self.grads['b'] = -1 / N * torch.matmul(torch.ones(size=[N]), (labels-self.outputs))

3.2.4.2 参数更新

在计算参数的梯度之后，我们使用3.1.4.2中实现的梯度下降法进行参数更新。

3.2.5 模型训练

实例化RunnerV2类，并传入训练配置。使用训练集和验证集进行模型训练，共训练500个epoch。每隔50个epoch打印训练集上的指标。

代码实现如下：

# 固定随机种子，保持每次运行结果一致
torch.manual_seed(102)

# 特征维度
input_dim = 2
# 类别数
output_dim = 3
# 学习率
lr = 0.1

# 实例化模型
model = model_SR(input_dim=input_dim, output_dim=output_dim)
# 指定优化器
optimizer = SimpleBatchGD(init_lr=lr, model=model)
# 指定损失函数
loss_fn = MultiCrossEntropyLoss()
# 指定评价方式
metric = accuracy
# 实例化RunnerV2类
runner = RunnerV2(model, optimizer, metric, loss_fn)

# 模型训练
runner.train([X_train, y_train], [X_dev, y_dev], num_epochs=500, log_eopchs=50, eval_epochs=1, save_path="best_model.pdparams")

# 可视化观察训练集与验证集的准确率变化情况
plot(runner,fig_name='linear-acc2.pdf')

运行结果：

best accuracy performence has been updated: 0.00000 --> 0.70625
[Train] epoch: 0, loss: 1.0986149311065674, score: 0.3218750059604645
[Dev] epoch: 0, loss: 1.0805636644363403, score: 0.706250011920929
best accuracy performence has been updated: 0.70625 --> 0.71250
best accuracy performence has been updated: 0.71250 --> 0.71875
best accuracy performence has been updated: 0.71875 --> 0.72500
best accuracy performence has been updated: 0.72500 --> 0.73125
best accuracy performence has been updated: 0.73125 --> 0.73750
best accuracy performence has been updated: 0.73750 --> 0.74375
best accuracy performence has been updated: 0.74375 --> 0.75000
best accuracy performence has been updated: 0.75000 --> 0.75625
best accuracy performence has been updated: 0.78750 --> 0.79375
best accuracy performence has been updated: 0.79375 --> 0.80000
[Train] epoch: 200, loss: 0.6921818852424622, score: 0.784375011920929
[Dev] epoch: 200, loss: 0.8020225763320923, score: 0.793749988079071
best accuracy performence has been updated: 0.80000 --> 0.80625
[Train] epoch: 300, loss: 0.6840379238128662, score: 0.7906249761581421
[Dev] epoch: 300, loss: 0.81141597032547, score: 0.8062499761581421
best accuracy performence has been updated: 0.80625 --> 0.81250
[Train] epoch: 400, loss: 0.680213987827301, score: 0.807812511920929
[Dev] epoch: 400, loss: 0.819807231426239, score: 0.8062499761581421

这里说一点，引用X包里的X函数，需要使用X.X，包名与包内函数重名的时候

3.2.6 模型评价

使用测试集对训练完成后的最终模型进行评价，观察模型在测试集上的准确率。

代码实现如下：

score, loss = runner.evaluate([X_test, y_test])
print("[Test] score/loss: {:.4f}/{:.4f}".format(score, loss))

运行结果：

[Test] score/loss: 0.8400/0.7014

可视化观察类别划分结果。

# 均匀生成40000个数据点
x1, x2 = torch.meshgrid(torch.linspace(-3.5, 2, 200), torch.linspace(-4.5, 3.5, 200), indexing='xy')
x = torch.stack([torch.flatten(x1), torch.flatten(x2)], axis=1)
# 预测对应类别
y = runner.predict(x)
y = torch.argmax(y, axis=1)
# 绘制类别区域
plt.ylabel('x2')
plt.xlabel('x1')
plt.scatter(x[:,0].tolist(), x[:,1].tolist(), c=y.tolist(), cmap=plt.cm.Spectral, label='before plot.scatter')

n_samples = 1000
X, y = make_multiclass_classification(n_samples=n_samples, n_features=2, n_classes=3, noise=0.2)

plt.scatter(X[:, 0].tolist(), X[:, 1].tolist(), marker='*', c=y.tolist(), label='after plot.scatter')
plt.legend()
plt.show()

运行结果：