机器学习练习2 python复现- 逻辑回归

在此练习中，需要实现逻辑回归应用于分类任务。还通过将正则化加入训练算法中来提高算法的鲁棒性，并用更复杂的情形进行测试。

逻辑回归

在训练的初始阶段，我们将要构建一个逻辑回归模型来预测，某个学生是否被大学录取。设想你是大学相关部分的管理者，想通过申请学生两次测试的评分，来决定他们是否被录取。现在你拥有之前申请学生的可以用于训练逻辑回归的训练样本集。对于每一个训练样本，你有他们两次测试的评分和最后是被录取的结果。为了完成这个预测任务，我们准备构建一个可以基于两次测试评分来评估录取可能性的分类模型。

数据检查

#导包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#数据导入与查看
path = 'ex2data1.txt'
data = pd.read_csv(path,header = None,names = ['Exam 1','Exam 2','Admitted'])
data.head()

#创建两次测试分数对应散点图，并使用颜色编码进行可视化，样本为正表示被录取，负表示未被录取
#取出被录取的两次测试数据
positive = data[data['Admitted'].isin([1])]
#取出未被录取的两次测试数据
negative = data[data['Admitted'].isin([0])]
#绘制散点图
fig,ax = plt.subplots(figsize = (12,8))
ax.scatter(positive['Exam 1'],positive['Exam 2'],s = 50,c = 'b',marker = 'o',label = 'Admitted')
ax.scatter(negative['Exam 1'],negative['Exam 2'],s = 50,c = 'r',marker = 'x',label = 'Not Admitted')
ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')
plt.show()

?

观察可以发现两类间有清晰的决策边界，现在我们需要实现逻辑回归，就可以训练一个模型来预测结果。

sigmoid函数

g代表一个常用的逻辑函数为S形函数，公式为：

合起来，我们得到逻辑回归模型的假设函数：

def sigmoid(z):
    return 1 / (1+np.exp(-z))

#sigmoid函数的简单运用
nums = np.arange(-10,10,step = 1)
fig,ax = plt.subplots(figsize = (12,8))
ax.plot(nums,sigmoid(nums),'r')
plt.show()

?

编写代价函数评估结果。

def cost(theta,X,y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y,np.log(sigmoid(X*theta.T)))
    second = np.multiply((1-y),np.log(1-sigmoid(X*theta.T)))
    return np.sum(first-second)/ (len(X))

#调整原始数据矩阵
data.insert(0,'Ones',1)

cols = data.shape[1]
X = data.iloc[:,0:cols - 1]
y = data.iloc[:,cols-1:cols]

X = np.array(X.values)
y = np.array(y.values)
theta = np.zeros(3)

theta

array([0., 0., 0.])

X.shape,theta.shape,y.shape

((100, 3), (3,), (100, 1))

计算初始化参数的代价函数(theta为0)

cost(theta,X,y)

0.6931471805599453

梯度下降(gradient descent)

def gradient(theta,X,y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X*theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error,X[:,i])
        grad[i] = np.sum(term) / len(X)
        
    return grad

注意，我们实际上没有在这个函数中执行梯度下降，我们仅仅在计算一个梯度步长。在练习中，一个称为“fminunc”的Octave函数是用来优化函数来计算成本和梯度参数。由于我们使用Python，我们可以用SciPy的“optimize”命名空间来做同样的事情。

使用数据和初始参数为0的梯度下降的结果。

gradient(theta,X,y)

array([ -0.1       , -12.00921659, -11.26284221])

使用SciPy's truncated newton（TNC）实现寻找最优参数。

import scipy.optimize as opt

result = opt.fmin_tnc(func = cost,x0=theta,fprime=gradient,args=(X,y))
result

(array([-25.16131872,   0.20623159,   0.20147149]), 36, 0)

代价函数计算的结果

cost(result[0],X,y)

0.20349770158947425

接下来，我们需要编写一个函数，用我们所学的参数theta来为数据集X输出预测。然后，我们可以使用这个函数来给我们的分类器的训练精度打分。 逻辑回归模型的假设函数：

当 ??? 大于等于0.5时，预测 y=1

当 ??? 小于0.5时，预测 y=0 。

def predict(theta,X):
    probability = sigmoid(X*theta.T)
    return [1 if x>= 0.5 else 0 for x in probability]

theta_min = np.matrix(result[0])
predictions = predict(theta_min,X)
correct = [1 if ((a==1 and b == 1) or (a==0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
accuracy = (sum(correct) % len(correct))
print('accuracy = {0}%'.format(accuracy))

accuracy = 89%

我们的逻辑回归分类器预测正确，如果一个学生被录取或没有录取，达到89%的精确度。不坏！记住，这是训练集的准确性。我们没有保持住了设置或使用交叉验证得到的真实逼近，所以这个数字有可能高于其真实值（这个话题将在以后说明）。

正则化逻辑回归

加入正则项提升逻辑回归算法。

设想你是工厂的生产主管，你有一些芯片在两次测试中的测试结果。对于这两次测试，你想决定是否芯片要被接受或抛弃。为了帮助你做出艰难的决定，你拥有过去芯片的测试数据集，从其中你可以构建一个逻辑回归模型。

path = 'ex2data2.txt'
data2 = pd.read_csv(path,header=None,names=['Test 1', 'Test 2', 'Accepted'])
data2.head()

positive = data2[data2['Accepted'].isin([1])]
negative = data2[data2['Accepted'].isin([0])]

fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['Test 1'], positive['Test 2'], s=50, c='b', marker='o', label='Accepted')
ax.scatter(negative['Test 1'], negative['Test 2'], s=50, c='r', marker='x', label='Rejected')
ax.legend()
ax.set_xlabel('Test 1 Score')
ax.set_ylabel('Test 2 Score')
plt.show()

?

注意到其中没有线性决策界限，来良好的分开两类数据。一个方法是用像逻辑回归这样的线性技术来构造从原始特征的多项式中得到的特征。让我们通过创建一组多项式特征入手吧。

degree = 5
x1 = data2['Test 1']
x2 = data2['Test 2']

data2.insert(3, 'Ones', 1)

for i in range(1, degree):
    for j in range(0, i):
        data2['F' + str(i) + str(j)] = np.power(x1, i-j) * np.power(x2, j)

data2.drop('Test 1', axis=1, inplace=True)
data2.drop('Test 2', axis=1, inplace=True)

data2.head()

此时需要修改第1部分的成本和梯度函数，包括正则化项，首先是代价函数

def costReg(theta,X,y,learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
    return np.sum(first - second) / len(X) + reg

添加正则化梯度函数：

如果我们需要梯度下降法令这个代价函数最小化，因为我们未对??0进行正则化，所以梯度下降分为两种情形：

整理可得：

def gradientReg(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        
        if (i == 0):
            grad[i] = np.sum(term) / len(X)
        else:
            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])
    
    return grad

初始化变量

# set X and y (remember from above that we moved the label to column 0)
cols = data2.shape[1]
X2 = data2.iloc[:,1:cols]
y2 = data2.iloc[:,0:1]

# convert to numpy arrays and initalize the parameter array theta
X2 = np.array(X2.values)
y2 = np.array(y2.values)
theta2 = np.zeros(11)

learningRate = 1

costReg(theta2, X2, y2, learningRate)

0.6931471805599454

gradientReg(theta2, X2, y2, learningRate)

array([0.00847458, 0.01878809, 0.05034464, 0.01150133, 0.01835599,
       0.00732393, 0.00819244, 0.03934862, 0.00223924, 0.01286005,
       0.00309594])

使用优化函数计算优化结果

result2 = opt.fmin_tnc(func=costReg, x0=theta2, fprime=gradientReg, args=(X2, y2, learningRate))
result2

(array([ 0.53010248,  0.29075567, -1.60725764, -0.58213819,  0.01781027,
        -0.21329508, -0.40024142, -1.37144139,  0.02264304, -0.9503358 ,
         0.0344085 ]),
 22,
 1)

最后使用第1部分中的预测函数查看我们方案在训练数据上的准确度

theta_min = np.matrix(result2[0])
predictions = predict(theta_min, X2)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y2)]
accuracy = (sum(map(int, correct)) % len(correct))
print ('accuracy = {0}%'.format(accuracy))

accuracy = 78%