吴恩达机器学习作业Python实现(二)：logistic回归

在这部分练习中，你将建立一个logistics回归模型来预测一个学生是否能被大学录取。假如你是大学招生办的工作人员，你想通过学生的两次考试成绩来决定他被录取的概率。你有一些往届学生的历史数据作为逻辑回归的训练集，对于每一个训练样本，你有学生两次考试的分数和录取结果。你的任务是建立一个分类模型来估计每个学生基于两次考试成绩被录取的概率。

Logistics回归

导入数据并查看

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.optimize as opt

path = r'ex2data1.txt'
data = pd.read_csv(path,names=['exam1','exam2','admitted'])
print(data.head())

       exam1      exam2  admitted
0  34.623660  78.024693         0
1  30.286711  43.894998         0
2  35.847409  72.902198         0
3  60.182599  86.308552         1
4  79.032736  75.344376         1

绘制散点图，横坐标为exam1的分数，纵坐标为exam2的分数，被录取则样本为正，没有被录取则样本为负，用颜色样式加以区分。
pandas库中的isin()函数用于数据筛选，isin()接受一个列表，判断该列中元素是否在列表中，多用于要选择某列等于多个数值或者字符串时。data[data[‘admitted’].isin([‘1’])]选取admitted列值为1的所有行，等价于data[data[‘admitted’]==1]。

positive = data[data.admitted.isin(['1'])]  # 1
negetive = data[data.admitted.isin(['0'])]  # 0

fig, ax = plt.subplots(figsize=(6,5))
ax.scatter(positive['exam1'], positive['exam2'], c='b', label='Admitted')
ax.scatter(negetive['exam1'], negetive['exam2'], s=50, c='r', marker='x', label='Unadmitted')
ax.legend()#默认将图例调整到最佳区域
# 设置横纵坐标名
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')
plt.show()

实现代价函数,变量初始化

logistic回归的假设函数为$$h_{\theta }=g({\theta }^{T}X)$$其中g是一个S型函数(Sigmoid function)$$g(z)=\frac{1}{1+e^{-z}}$$
逻辑回归的代价函数如下，这个函数是凸的。$$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)}))]$$
此过程应用向量化方法，矩阵相乘时需要考虑矩阵的维度，必要时要对矩阵进行转置。注意矩阵matrix和数组array的区别。

def Sigmoid(z):
    return 1/(1+np.exp(-z))
def computeCost(theta, X, Y):
    first = (-Y) * np.log(Sigmoid(X @ theta))
    second = (1 - Y)*np.log(1 - Sigmoid(X @ theta))
    return np.mean(first - second)

# add a ones column - this makes the matrix multiplication work out easier
if 'Ones' not in data.columns:
    data.insert(0, 'Ones', 1)
    
# set X (training data) and y (target variable)
X = data.iloc[:, :-1].values  # Convert the frame to its Numpy-array representation.
Y = data.iloc[:, -1].values  # Return is NOT a Numpy-matrix, rather, a Numpy-array.
theta = np.zeros(X.shape[1])
print(computeCost(theta, X, Y))

0.6931471805599453

计算梯度

bathch gradient decent(批量梯度下降)
$${\theta }_j:={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}J(\theta )$$
$${\theta }_j:={\theta }_j-\alpha \frac{{\sum }_{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}){x_j}^{(i)}}{m}$$
实现梯度下降过程中(如下)程序报错ValueError: shapes (100,1) and (3,100) not aligned: 1 (dim 1) != 3 (dim 0),实在不知道哪里出错了，多次检查矩阵乘法未果，遂放弃梯度下降算法。在此只计算出梯度：$$\frac{\partial }{\partial {\theta }_j}J(\theta )$$

"""
def gradientDescent(X, Y, theta, alpha, epoch):
    temp = np.matrix(np.zeros(theta.shape))  # 初始化一个临时矩阵(1, 3)
    cost = np.zeros(epoch)                   # 初始化一个ndarray，包含每次epoch的cost
    m = X.shape[0]                           # m为总的样本数
    #利用向量化一步求解
    for i in range(epoch):
        temp = theta - (alpha / m) * (Sigmoid(X * theta.T)-Y).T * X
        theta = temp
        cost[i] = computeCost(X,Y,theta) #得到每次迭代代价函数的值
    return theta,cost

#初始化学习率α和要进行迭代的次数
alpha = 0.01
epoch = 1000
#现在让我们运行梯度下降算法来将我们的参数θ适合于训练集
final_theta,cost = gradientDescent(X,Y,theta,alpha,epoch)
print(cost)
"""
def gradient(theta, X, y):
    return (X.T @ (Sigmoid(X @ theta) - y))/len(X)   
# the gradient of the cost is a vector of the same length as θ where the jth element (for j = 0, 1, . . . , n)
print(gradient(theta, X, Y))

[ -0.1        -12.00921659 -11.26284221]

运用Scipy优化算法来最小化代价函数，学习参数$\theta$

由于上面并没有实现梯度下降算法，只是计算出梯度。在练习中，一个称为“fminunc”的Octave函数是用来优化函数来计算成本和梯度参数。由于我们使用Python，我们可以用SciPy的“optimize”命名空间来做同样的事情。scipy中的optimize子包中提供了常用的最优化算法函数实现，我们可以直接调用这些函数完成我们的优化问题。在用python实现逻辑回归和线性回归时，使用梯度下降法最小化cost function，用到了fmin_tnc()和minimize()。具体用法可以参考https://www.cnblogs.com/tongtong123/p/10634716.html
这里我们使用的是高级优化算法，运行速度通常远远超过梯度下降。方便快捷。
只需传入cost函数，已经所求的变量theta，和梯度。cost函数定义变量时变量tehta要放在第一个，若cost函数只返回cost，则设置fprime=gradient。

1.fmin_tnc()

scipy.optimize.fmin_tnc(func, x0, fprime=None, args=(), approx_grad=0)，这里只列举出常用参数，详细可参考官方文档https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.optimize.fmin_tnc.html
参数：
func：优化的目标函数（此例中为代价函数）
x0(array_like)：初值（此例中为$\theta$）
fprime：提供优化函数func的梯度函数，不然优化函数func必须返回函数值和梯度，或者设置approx_grad=True
approx_grad :如果设置为True，会给出近似梯度
args：元组，是传递给优化函数的参数
返回值：
x(ndarray): The solution.
nfeval(int): The number of function evaluations.
rc(int): Return code, see below

import scipy.optimize as opt
result = opt.fmin_tnc(func=computeCost, x0=theta, fprime=gradient, args=(X, Y))
print(result[0])   #返回代价函数取最小值时的参数值

[-25.1613186    0.20623159   0.20147149]

2.minimize()

除了用fmin_tnc()方法，还可以用minimize()方法，
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None)，这里只列举出常用参数，详细可参考官方文档https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
参数：
fun ：优化的目标函数
x0 ：初值，一维数组，shape (n,)
args ： 元组，可选，额外传递给优化函数的参数
method：求解的算法，选择TNC则和fmin_tnc()类似
jac：返回梯度向量的函数
返回值：
返回优化结果对象，x：优化问题的目标数组。success: True表示成功与否，不成功会给出失败信息。

res = opt.minimize(fun=computeCost, x0=theta, args=(X, Y), method='TNC', jac=gradient)
print(res.x)

[-25.1613186    0.20623159   0.20147149]

输出预测值及被录取（输出值为1）的概率

逻辑回归模型的假设函数：$$h_{\theta }(x)=\frac{1}{1+e^{{\theta }^{T}X}}$$
$h_{\theta }(x)$的作用是，对于给定的输入变量，根据选择的参数计算输出变量=1 的可能性。
当$h_{\theta }(x)$小于0.5时，预测 y = 0
当$h_{\theta }(x)$大于0.5时，预测 y = 1
接下来，我们需要编写一个函数，用我们所学的参数theta来为数据集X输出预测。然后，我们可以使用这个函数来给我们的分类器的训练精度打分。

#对样本集进行预测，得到一个预测值列表
def predict(theta, X):
    probability = Sigmoid(X@theta)
    return [1 if x >= 0.5 else 0 for x in probability]  # return a list

final_theta = result[0]
predictions = predict(final_theta, X)
correct = [1 if a==b else 0 for (a, b) in zip(predictions, Y)]

#计算预测精度
accuracy = sum(correct) / len(X)
print(accuracy)

0.89

绘制决策边界

分隔 y=0 的区域和 y=1 的区域的曲线即为决策边界，在这个假设模型中是${\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2=0$所对应的曲线。

x1 = np.arange(130, step=0.1)
x2 = -(final_theta[0] + x1*final_theta[1]) / final_theta[2]

fig, ax = plt.subplots(figsize=(8,5))
ax.scatter(positive['exam1'], positive['exam2'], c='b', label='Admitted')
ax.scatter(negetive['exam1'], negetive['exam2'], s=50, c='r', marker='x', label='Not Admitted')
ax.plot(x1, x2)
ax.set_xlim(0, 130)
ax.set_ylim(0, 130)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_title('Decision Boundary')
plt.show()

正则化Logistic回归

在训练的第二部分，我们将要通过加入正则项提升逻辑回归算法。简而言之，正则化是成本函数中的一个术语，它使算法更倾向于“更简单”的模型（在这种情况下，模型将更小的系数）。这个理论助于减少过拟合，提高模型的泛化能力。
在练习的这一部分中，你将实现正则化的逻辑回归，以预测来自制造工厂的微芯片是否通过质量保证(QA)。在这过程中，每个芯片都要经过各种测试，以确保其正常工作。假设你是工厂的产品经理，你有一些微芯片在两种不同的测试上的测试结果。从这两个测试中，你想确定应该接受还是拒绝微芯片。为了帮助你做出决策，你有一个关于过去微芯片的测试结果的数据集，你可以从中构建一个逻辑回归模型。

导入数据并查看

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

	Test 1	Test 2	Accepted
0	0.051267	0.69956	1
1	-0.092742	0.68494	1
2	-0.213710	0.69225	1
3	-0.375000	0.50219	1
4	-0.513250	0.46564	1

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

fig, ax = plt.subplots(figsize=(8,5))
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()

注意到正负样本之间没有线性的决策边界，所以直接用Logistic回归在这个数据集上并不能表现良好。

生成多项式特征

一个拟合数据的更好的方法是从每个数据点创建更多的特征。
我们将把这些特征映射到所有的x1和x2的多项式项上，直到第六次幂。

def feature_mapping(x1, x2, power):
    data = {}
    for i in np.arange(power + 1):
        for p in np.arange(i + 1):
            data["f{}{}".format(i - p, p)] = np.power(x1, i - p) * np.power(x2, p)
    return pd.DataFrame(data)
x1 = data2['Test 1'].values
x2 = data2['Test 2'].values
_data2 = feature_mapping(x1, x2, power=6)
_data2.head()

	f00	f10	f01	f20	f11	f02	f30	f21	f12	f03	...	f23	f14	f05	f60	f51	f42	f33	f24	f15	f06
0	1.0	0.051267	0.69956	0.002628	0.035864	0.489384	0.000135	0.001839	0.025089	0.342354	...	0.000900	0.012278	0.167542	1.815630e-08	2.477505e-07	0.000003	0.000046	0.000629	0.008589	0.117206
1	1.0	-0.092742	0.68494	0.008601	-0.063523	0.469143	-0.000798	0.005891	-0.043509	0.321335	...	0.002764	-0.020412	0.150752	6.362953e-07	-4.699318e-06	0.000035	-0.000256	0.001893	-0.013981	0.103256
2	1.0	-0.213710	0.69225	0.045672	-0.147941	0.479210	-0.009761	0.031616	-0.102412	0.331733	...	0.015151	-0.049077	0.158970	9.526844e-05	-3.085938e-04	0.001000	-0.003238	0.010488	-0.033973	0.110047
3	1.0	-0.375000	0.50219	0.140625	-0.188321	0.252195	-0.052734	0.070620	-0.094573	0.126650	...	0.017810	-0.023851	0.031940	2.780914e-03	-3.724126e-03	0.004987	-0.006679	0.008944	-0.011978	0.016040
4	1.0	-0.513250	0.46564	0.263426	-0.238990	0.216821	-0.135203	0.122661	-0.111283	0.100960	...	0.026596	-0.024128	0.021890	1.827990e-02	-1.658422e-02	0.015046	-0.013650	0.012384	-0.011235	0.010193

5 rows × 28 columns

经过映射，我们将有两个特征的向量转化成了一个28维的向量。
在这个高维特征向量上训练的logistic回归分类器将会有一个更复杂的决策边界，当我们在二维图中绘制时，会出现非线性。
虽然特征映射允许我们构建一个更有表现力的分类器，但它也更容易过拟合。在接下来的练习中，我们将实现正则化的logistic回归来拟合数据，并且可以看到正则化如何帮助解决过拟合的问题。

实现正则化代价函数，变量初始化

$$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)}))]$$
注意：不惩罚第一项${\theta }_0$

def costReg(theta, X, Y, lambuda=1):
    # 不惩罚第一项
    _theta = theta[1: ]
    reg = (lambuda / (2 * len(X))) *(_theta @ _theta)  # _theta@_theta == inner product
    return computeCost(theta, X, Y) + reg

# 这里因为做特征映射的时候已经添加了偏置项，所以不用手动添加了。
X = _data2.values  
Y = data2['Accepted'].values
theta = np.zeros(X.shape[1])
print(costReg(theta,X,Y,1))

0.6931471805599454

正则化梯度

正则化梯度下降：
$${\theta }_0:={\theta }_0-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}){x_0}^{(i)}$$
$${\theta }_j:={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m[(h_{\theta }(x^{(i)})-y^{(i)}){x_j}^{(i)}+\frac{\lambda }{m}{\theta }_j],forj=1,2,...n$$
对上面的算法中 j=1,2,…,n 时的更新式子进行调整可得：
$${\theta }_j:={\theta }_j(1-a\frac{\lambda }{m})-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}){x_j}^{(i)}$$
仅需要计算出梯度：$$\frac{\partial }{\partial {\theta }_j}J(\theta )$$

def gradientReg(theta, X, Y, lambuda=1):
    reg = (lambuda / len(X)) * theta
    reg[0] = 0  
    return gradient(theta, X, Y) + reg

print(gradientReg(theta,X,Y,1))

[8.47457627e-03 1.87880932e-02 7.77711864e-05 5.03446395e-02
 1.15013308e-02 3.76648474e-02 1.83559872e-02 7.32393391e-03
 8.19244468e-03 2.34764889e-02 3.93486234e-02 2.23923907e-03
 1.28600503e-02 3.09593720e-03 3.93028171e-02 1.99707467e-02
 4.32983232e-03 3.38643902e-03 5.83822078e-03 4.47629067e-03
 3.10079849e-02 3.10312442e-02 1.09740238e-03 6.31570797e-03
 4.08503006e-04 7.26504316e-03 1.37646175e-03 3.87936363e-02]

运用Scipy优化算法来最小化代价函数，学习参数

result2 = opt.fmin_tnc(func=costReg, x0=theta, fprime=gradientReg, args=(X, Y, 2))
print(result2[0])

[ 0.90267454  0.33721089  0.76006404 -1.39757946 -0.51417075 -0.91389985
  0.01516214 -0.21926017 -0.22677642 -0.16219637 -1.01270257 -0.04169398
 -0.39984069 -0.14458017 -0.82296284 -0.20346048 -0.13186937 -0.04837714
 -0.17183934 -0.17077936 -0.38820995 -0.72773035  0.00607685 -0.19391899
  0.00314606 -0.21203169 -0.06947222 -0.69320886]

输出预测值及被录取（输出值为1）的概率

final_theta = result2[0]
predictions = predict(final_theta, X)
correct = [1 if a==b else 0 for (a, b) in zip(predictions, Y)]
accuracy = sum(correct) / len(correct)
#输出预测精度
print(accuracy)

0.8305084745762712

高级Python库scikit-learn

from sklearn import linear_model#调用sklearn的线性回归包
model = linear_model.LogisticRegression(penalty='l2', C=1.0)
model.fit(X, Y.ravel())
print(model.score(X,Y))

0.8305084745762712

绘制决策边界

x = np.linspace(-1, 1.5, 250)
xx, yy = np.meshgrid(x, x)

z = feature_mapping(xx.ravel(), yy.ravel(), 6).values
z = z @ final_theta
z = z.reshape(xx.shape)

plt.scatter(positive['Test 1'], positive['Test 2'], s=50, c='b', marker='o', label='Accepted')
plt.scatter(negative['Test 1'], negative['Test 2'], s=50, c='r', marker='x', label='Rejected')
plt.contour(xx, yy, z, 0)
plt.ylim(-.8, 1.2)

(-0.8, 1.2)