逻辑回归在多分类问题中的应用

在前面我写过用KNN算法处理数字的手写识别，但是KNN算法是“懒惰学习”的代表，不具有显式学习过程。此类学习技术在训练过程中仅仅是把样本保存起来，训练时间为0，待收到训练样本后在进行处理。但是这类算法不产生固定的参数，往往对于每一个新样本都需要较长计算过程。所以这次使用逻辑回归处理此类问题。

首先导入需要使用的模块，导入loadmat读取mat类型文件，以及numpy对数据进行处理：

import numpy as np
from scipy.io import loadmat

然后我们查看并处理数据：

data = loadmat('F:\\MachineLearning\data\ex3data1.mat')
print(data)

我们可以看到一个字典，其中包含样本信息的矩阵X和样本标记二维数组y：

{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',
 '__version__': '1.0',
 '__globals__': [],
 'X': array([[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]),
 'y': array([[10],
        [10],
        [10],
        ...,
        [ 9],
        [ 9],
        [ 9]], dtype=uint8)}

所以我问可以分别取X、y，并查看其数据类型：

X = data['X']
m,d = X.shape
y = data['y']
y = y.reshape((m,))
X = np.insert(X,0,values=np.ones(m),axis=1)
print(X.shape,y.shape)

输出如下所示：

(5000, 401) (5000,)

即我们的训练集中共有5000个样本，每个样本的维度为401，由于我们在回归过程中需要处理常数项，所以在X矩阵中第一列插入全为1的数组。除去第一列，可以知道每个样本具有400个属性。第四行中，我们把y从二维数组形式转化为一维数组，方便我们下面的处理。

下面我们遍历y中的每个类别：

y_set = set(list(y))
for i in y_set:
	print(i)

得到：

1
2
3
4
5
6
7
8
9
10

即样本类别为1到10。

下面为逻辑回归中的sigmoid函数,代价函数及梯度函数。这里为了防止过拟合，采用正则化形式的代价和梯度函数：
sigmoid函数$$g(z)=\frac{1}{1+e^{-z}}$$ $$h_{\theta }(x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$
正则化逻辑回归代价函数$$J(\theta )=\frac{1}{m}\sum _{i=1}^m[-y_i\log (h_{\theta }(x_i))-(1-y_i)\log (1-h_{\theta }(x_i)]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$
正则化逻辑回归梯度函数
对于$i\ge 1$
$$\frac{\partial J(\theta )}{\partial {\theta }_j}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x_i)-y_i)x_{ij}+\frac{\lambda }{m}{\theta }_j$$
对于$i=0$
$$\frac{\partial J(\theta )}{\partial {\theta }_0}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x_i)-y_i)x_{i0}$$

所以下面我们利用上面公式编写代价函数和梯度函数：

def regularized_cost(theta,X,y,regularized_param):
	m = X.shape[0]
	sigVector = sigmoid(X @ theta)
	cost_term = (-y.T @ sigVector - (1-y).T @ (1-sigVector))/m
	regularized_term = (regularized_param / (2*m)) * (tehta[1:].T @ theta[1:])
	cost = cost_term + regularized_term
	return cost

def regularized_gradient(theta,X,y,regularized_term):
	m = X.shape[0]
	sigVector = sigmoid(X @ theta)
	deltaVector = sigVector - y
	gradient_term = (X.T @ deltaVector)/m
	theta1 = np.insert(theta[1:],0,values=0)
	regularized_term = (regularized_term / m) * tehta1
	gradient = gradient_term + regularized_term
	return gradient

编写完代价函数和梯度函数，我们就可以利用scipy中minimize对模型参数进行估计，这里我们以此将1到10的类别作为正例，其他的作为反例，这样我们对每个样本都可以得到一组参数，我们挑选使此样本的sigmoid函数最大参数对应的类别最为预测的类别：

from scipy.optimize import minimize
def one_vs_all(X,y,regularized_param):
	m,d = X.shape
	num_label = len(set(list(y)))                           #类别数
	all_theta = np.zeros((num_label，d))
	for i in rang(1,num_label):
		theta=np.zeros(d)
		y_i = [1 if i==label else 0 for label in y]
		y_i = np.array(y_i)
		fmin = minimize(fun = regularized_cost,x0=theta,args=(X,y_i,regularized_term),method='TNC',\
			jac=regularized_gradient)
		all_theta[i-1] = fmin.x
	return all_theta

将X，y带入，并将正则化常数设为1：

print(one_vs_all(X,y,1))

得到将每一个类作为正例，其他作为反例的模型参数组成的矩阵，即每一行为每一类别作为正例的参数，一共有十个类别，所以有10行。得到：

[[-2.38187334e+00  0.00000000e+00  0.00000000e+00 ...  1.30433279e-03
  -7.29580949e-10  0.00000000e+00]
 [-3.18303389e+00  0.00000000e+00  0.00000000e+00 ...  4.46340729e-03
  -5.08870029e-04  0.00000000e+00]
 [-4.79638233e+00  0.00000000e+00  0.00000000e+00 ... -2.87468695e-05
  -2.47395863e-07  0.00000000e+00]
 ...
 [-7.98700752e+00  0.00000000e+00  0.00000000e+00 ... -8.94576566e-05
   7.21256372e-06  0.00000000e+00]
 [-4.57358931e+00  0.00000000e+00  0.00000000e+00 ... -1.33390955e-03
   9.96868542e-05  0.00000000e+00]
 [-5.40542751e+00  0.00000000e+00  0.00000000e+00 ... -1.16613537e-04
   7.88124085e-06  0.00000000e+00]]

下面对我们的训练集的的正确率进行预测：
对一个样本x,利用all_theta进行预测

def predict(all_theta,x):
	h = sigmoid(all_theta @ x)
	max_index =np.argmax(h)
	max_h = max_index + 1
	return max_h

计算训练集正确率

def predict_all(X,all_theta,y):
	num_correct=0
	for i in range(m):
		x=X[i,:]
		predict_label=predict(all_theta,x)
		if predict_label == y[i]:
			num_correct+=1
	accuracy = num_correct/m
	return accuracy

下面计算正确率：

print(predict_all(X_all_theta,y))

得到：

0.9446