spss多元无序logistic回归模型word

 

本数据集为.mat文件，包含两个变量 $X,y$

变量$X$包含了5000个手写数字的样本图片，每个样本图片大小为20×20像素，因此将每个样本表示为400维，其中的数值代表每个像素格的灰度值。

变量$y$则是5000个样本所对应的标签，例如1-500的样本$y = 10$，表示该500个样本是对应0类，501-1000的样本$y = 1$，表示该500个样本是对应1类，如此下去，$y$的取值从1-10，分别对应数字1-9和0

1.多元分类

 

所谓的多元分类就是，可以理解为：二分类的再分类问题。

例如：我们想要将数据集A分类称为a，b，c，d四类。首先我们可以将A分成a和B两类，之后将B分成b和C两类，最后将C分类为c和d两类。

这种多元分类问题是基于二分类实现的，可以称为one vs others分类。

多元分类问题，其实就是使用多个二分类器实现的，其在本质上同逻辑回归二分类问题一致。

2.标签向量化

如上图所示：我们从原始的数据集中获得的标签数据是一个（5000,1）的列向量。

同样，我们参考二分类问题。二分类问题的标签仅存在0和1两类。

即将A分成a和B两类时，我们可以令a所对应的标签转化为1，其他数据对应标签转换为0，成功分离出a。此时获取到一个新的标签向量$y_1$，仅含有0和1。

同理，我们可以使用同样的办法，分离出b，c，d。此时则应该有三个新的标签向量$y_2,y_3,y_4$。

这个过程就称为标签的向量化过程。

3.多元分类流程

• 读取数据
• 获取特征值$X$
• 获取标签值$y$
• 标签值向量化$y[0],y[1],y[2],y[3]...$
• 特征缩放(如有必要)
• 构建特征与标签之间的函数关系（假设函数）
• 构建正则化代价函数
• 求取正则化代价函数的偏导数
• 最小正则化化代价函数，并获取到此时的权重
• 根据权重得到特征与标签之间的函数关系

多元分类的流程相比于二分类仅仅多了一个标签值向量化，其余和二分类完全相同。需要注意的是，这里的$y[0]$就相当于一个二分类器，其对应一组$θ[0]$，而$y[1]$又相当于另一个二分类器，其又对应一组$θ[1]$，如此，每一个标签向量都相当于一个二分类器，对应一组权重系数。

4.(n,) (n,1)之间的区别

在进行代码实现之前需要明确一个问题。这将有助于我们理解代码中涉及计算相关的维度对应问题。

numpy中，关于(n,) (n,1)之间的区别

我们可以观察以下代码：

import numpy as np

list_1 = [[1, 2, 3]]
a = np.array(list_1)
print(a)  # [[1 2 3]]
print(a.shape)  # (1, 3)
b = np.arange(1, 4, 1).reshape(3, 1)
print(b)
# [[1]
#  [2]
#  [3]]
print(b.shape)   #(3, 1)
print(a @ b)  # [[14]]
# 这里的计算是(1,3)(3,1)维度对应的计算，这很符合我们的矩阵计算规则

import numpy as np

list_1 = [[1, 2, 3]]
a = np.array(list_1).reshape(3)
print(a)  # [1 2 3]
print(a.shape)  # (3,)
b = np.arange(1, 7, 1).reshape(3, 2)
print(b)
# [[1]
#  [2]
#  [3]]
print(b.shape)  # (3, 1)
print(a @ b)   # [14]
# 我们发现，这段代码给出了与上段代码相似的结果，这是(3,)(3,1)维度对应的计算，如果基于矩阵计算规则，这是否可以让我们有如下的理解：
# (3,)等价于(1,3)的行向量？

import numpy as np

list_1 = [[1, 2, 3]]
a = np.array(list_1).reshape(3)
print(a)  # [1 2 3]
print(a.shape)  # (3,)
b = np.arange(1, 4, 1).reshape(1,3)
print(b)  # [[1 2 3]]
print(b.shape)  # (1, 3)
print(b @ a)   # [14]
# 在这段代码中，我们得到了与第二段代码相同的结果，这是(1,3)(3，)维度对应的计算，如果基于矩阵计算规则，这是否可以让我们有如下的理解：
# (3,)等价于(3,1)的列向量？
# 我们试图计算(3,)(1,3)的时候，发现存在语法错误

通过上述代码，我们获得了如下的结论：

**(n,)在进行矩阵运算的时候，其可以根据情况选择表现为列向量还是行向量。而这个情况就是，其在遵循矩阵计算规则的前提下，计算结果必须是一维数组，否则将会报错。**因此将(n,1)转化成(n,)可以在一定程度上忽略维度对应。

5.代码实现

5.1数据图形化

变量$X$包含了5000个手写数字的样本图片，每个样本图片大小为20×20像素，因此将每个样本表示为400维，其中的数值代表每个像素格的灰度值。

变量$y$则是5000个样本所对应的标签

为了更好的理解数据的真实含义，可以通过如下代码将数据图形化：

import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
import matplotlib


def load_data(path):  # 加载数据
    data = sio.loadmat(path)  # 通过路径读取.mat数据
    y = data.get('y')  # 获取到标签y(5000,1)
    y = y.reshape(y.shape[0])  # 为了方便计算，将y转换成(5000,)的一维数组
    X = data.get('X')  # (5000,400)
    return X, y


def plot_an_image(image):  # 绘图
    fig, ax = plt.subplots(figsize=(1, 1))  # 定义一张大小为(1,1)的图像
    ax.matshow(image.reshape((20, 20)), cmap=matplotlib.cm.binary)  # 把矩阵或数组绘制成图像的函数
    plt.xticks(np.array([]))   # 不显示x轴
    plt.yticks(np.array([]))   # 不显示y轴


X, raw_y = load_data('ex3data1.mat')  # 读取特征X
pick_one = np.random.randint(0, 5000)  # 随机挑选一个样本
print(pick_one)
plot_an_image(X[pick_one, :])  # 根据挑选样本的特征绘图，即根据像素格的灰度值绘图
plt.show()
print('this should be {}'.format(raw_y[pick_one]))

运行结果：

4004   # 即第4004个样本
this should be 8

5.2完整代码

 

import numpy as np
import scipy.io as sio
import scipy.optimize as opt
from sklearn.metrics import classification_report  # 这个包是评价报告


def load_data(path):  # 加载数据
    data = sio.loadmat(path)  # 通过路径读取.mat数据
    y = data.get('y')  # 获取到标签y(5000,1)
    y = y.reshape(y.shape[0])  # 为了方便计算，将y转换成(5000,)的一维数组
    raw_X = data.get('X')  # (5000,400)
    X = np.insert(raw_X, 0, values=np.ones(raw_X.shape[0]), axis=1)  # 插入了第一列（全部为1），即x_0 = 1
    return X, y


def label_vectorization(raw_y):  # 标签向量化
    y_matrix = []  # 定义一个空列表，用来存储向量化后的数据
    for k in range(1, 11):  # 标签：1-10，对应1-9和0，将标签向量化为10个向量
        y_matrix.append((raw_y == k).astype(int))
    y_matrix = [y_matrix[-1]] + y_matrix[:-1]  # 由于最后一个向量对应的是0类，因此将0类提到第一个，即0-9
    y = np.array(y_matrix)  # 将list数据类型，转换成ndarray数据类型
    return y


def sigmoid(z):  # 定义sigmoid函数
    return 1 / (1 + np.exp(-z))


def cost(theta, X, y):  # 定义代价函数，这里需要注意X,theta,y的维度，其并不影响计算
    return np.mean(-y * np.log(sigmoid(X @ theta)) - (1 - y) * np.log(1 - sigmoid(X @ theta)))
    # 值得一提的是，在numpy中，+，-，*，/都只是相对应元素的运算，不遵循矩阵计算规则
    # numpy中的矩阵计算需要使用@或者.dot()函数
    # 同样也可以将ndarray类型数据通过np.matrix()转换成矩阵，使用*进行计算


def regularized_cost(theta, X, y, l=1):  # 正则化代价函数  l 为正则系数
    theta_j1_to_n = theta[1:]  # 从θ_1到θ_n ，通常情况下，我们不对θ_0进行正则化
    regularized_term = (l / (2 * len(X))) * np.power(theta_j1_to_n, 2).sum()  # 代价函数正则项，是一个数值
    return cost(theta, X, y) + regularized_term


def gradient(theta, X, y):  # 求偏导
    return (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y)


def regularized_gradient(theta, X, y, l=1):  # 正则化偏导  l 为正则系数
    theta_j1_to_n = theta[1:]  # # 从θ_1到θ_n
    regularized_theta = (l / len(X)) * theta_j1_to_n  # 代价函数偏导正则项，是一个数组
    # 由于没有正则化θ_0，因此拼接一个数组0，保证数组维数
    regularized_term = np.concatenate([np.array([0]), regularized_theta])
    return gradient(theta, X, y) + regularized_term


def logistic_regression(X, y, l=1):
    theta = np.zeros(X.shape[1])  # 定义θ 要和特征数据维度相同
    # 使用优化器中的优化函数，采用的方法是TNC“截断牛顿法”
    res = opt.minimize(fun=regularized_cost,
                       x0=theta,
                       args=(X, y, l),
                       method='TNC',
                       jac=regularized_gradient)
    final_theta = res.x  # 获取最终的θ值
    return final_theta


def predict(x, theta):  # 预测函数
    prob = sigmoid(x @ theta)
    return (prob >= 0.5).astype(int)


# 获取特征X，标签y
X, raw_y = load_data('ex3data1.mat')
y = label_vectorization(raw_y)  # 标签向量化

# 逻辑回归
k_theta = np.array([logistic_regression(X, y[k]) for k in range(10)])
# y[k].shape = (5000,)可视为行向量和列向量
# k_theta.shape = (10,401)

# 预测和验证
prob_matrix = sigmoid(X @ k_theta.T)  # h_θ（x），由于是多个二分类器组成的，故应是一个矩阵
# (5000,401)[(10,401).T]===>(5000,10),表示5000个样本，每个样本属于0-9的可能性，可能性最大的值，即为预测值。

y_pred = np.argmax(prob_matrix, axis=1)  # 返回沿轴axis最大值的索引，axis=1代表行.索引为0-9，更好对应分类

y_ture = raw_y.copy()  # 未进行向量化的标签1-9和10
y_ture[y_ture == 10] = 0  # 将标签中的10替换为0
print(classification_report(y_ture, y_pred))  # 评价报告

6.运行结果

precision    recall  f1-score   support

           0       0.97      0.99      0.98       500
           1       0.95      0.99      0.97       500
           2       0.95      0.92      0.93       500
           3       0.95      0.91      0.93       500
           4       0.95      0.95      0.95       500
           5       0.92      0.92      0.92       500
           6       0.97      0.98      0.97       500
           7       0.95      0.95      0.95       500
           8       0.93      0.92      0.92       500
           9       0.92      0.92      0.92       500

    accuracy                           0.94      5000
   macro avg       0.94      0.94      0.94      5000
weighted avg       0.94      0.94      0.94      5000

对0，6的识别正确率很高，达到98%，97%，而对5和9的识别率较低，只有92%。整体识别率达到94%。