专栏推荐

基础知识准备

先来看一下逻辑回归的假设函数：Z等于θTX

逻辑回归的代价函数为：

向量化的逻辑回归的损失函数（矩阵形式）：向量形式就是X（i），Y（i）都表示一个样本，而向量X和Y表示的是整个数据集X和标注Y

其中g（θX）就是假设函数h

有了代价函数之后，我们下面的任务是求解代价函数的最小值，多对应的参数θ，这里我们使用梯度下降的方式来完成这个操作。

偏导数的求导过程为：

首先将J（θ）分成以下几个部分，然后逐步传递求导，这就是复合求偏导

θTX=θ1X1+θ2X2+……+θmXm

对哪个θi求偏导，哪个θi就是变量，其它的θ就是常数，X也是常熟

其中（1/1+e-θTX）的倒数，可以看成是（a/b）的倒数形式

不要忘记θTX=θ1X1+θ2X2+……+θmXm，通过这个式子就可以求出对θj的偏导了

求偏导结果（梯度）：

向量化的偏导结果（梯度）：

我们完成梯度下降之后，就可以获取到θ,进而获取到假设函数，然后就可以通过假设函数来预测其它的数据集了。

假设函数代码为：

参数为Z，只要调用方法的时候，传递θTX，或者X和θ的点积就可以返回假设函数了。

代价函数代码为：

M为样本的数量，h为假设函数，dot为矩阵乘法，np.log为数学中的log，T为转秩

我们获取到代价函数J之后，在处理数据时遇到NAN值（缺失值）的几率还是比较大的，有的时候需要对数据值是否为nan值做判断，如果是的话，那么就返回np.inf,这个表示无限大的正数。我们的这个代价函数不是为了求偏导，因为求偏导我们直接用

这个公式就ok，所以也就是说我们可以直接跳过代价函数，那我们写代价函数的意义就是为了画代价函数的图像，通过可视化的方式来判断J的梯度下降的过程。

这是梯度下降（向量化）的求偏导数的方法，这个是同时对所有的θ求偏导，grad.flatten()方法为将grad折叠成的一维数组返回。注意1.0不是1，幻风就是在这里出现了问题。

逻辑回归的执行过程是这样的：我们不停的调用梯度下降函数，进行梯度下降，直到代价函数变化很小了，那此时的θ就是我们所要使用的θ

代码实现

我们现在有一个这样的数据集：data1.txt

它有三列，前两列为特征，第三列为标签。分为两类0和1

data = np.loadtxt('data1.txt', delimiter=',')

使用该代码可以实现读取文件data1.txt，然后使用,为分割符，输出data结果如下所示：整个数据集为100行3列

现在我们要对数据进行处理，分成样本X和标签y

X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]y = np.c_[data[:,2]]y = np.c_[data[:,2]]

在x上添加一列x0，为x0=1，np.c为行数一样的矩阵进行合并

data[0:2]是获取data数据的0列和第1列。

initial_theta = np.zeros(X.shape[1])

这个代码是初始化theta

cost = costFunction(initial_theta, X, y)grad = gradient(initial_theta, X, y)

这个代码的意义不大就是看初始theta的情况下此时的代价是多少，以及初始theta在梯度下降一次时的theta（grad）变成了什么。

res = minimize(costFunction, initial_theta, args=(X,y), method=None, jac=gradient, options={'maxiter':400})

这个代码的意思就是使用高级函数minimize来完成梯度下降，第一个参数是代价函数，第二个参数是初始化参数，第三个参数是样本数据，method有以下几个参数：CG，Newton-CG，BFGS，L-BFGS-B。参数maxiter的意思是最大梯度下降的次数。

我们并没有自己写循环来来回回的调用梯度下降，而是直接调用现成的minimize的方法来训练我们的模型，最终返回res是我们训练好的模型信息，而res.x的就是此时的已经训练好的最终的参数theta

那么此时我们就可以用这个参数来预测我们的其它数据了，我们先来写一个预测函数

def predict(theta, X, threshold=0.5): p = sigmoid(X.dot(theta.T)) >= threshold return(p.astype('int'))

这个测试函数是调用sigmoid函数，如果值大于0.5，那就表示为1，如果小于0.5那就等于0，下面预测

a=np.array([1, 10, 10])b=predict(res.x,a)

这个就是预测x1=1，x2=50，x3=50这个样本，那么在这个模型的情况下，b就是预测结果，值为：0，就表示此时类别为0，就表示在科目一为10分，科目二为10分的情况下成绩不合格。

我们如何计算出我们这个模型的准确度是多少呢？

方法就是我们将我们训练数据带入到训练模型中，然后看预测结果和实际结果一致的情况

p = predict(res.x, X)

print('Train accuracy {}%'.format(100*sum(p == y.ravel())/p.size))

这个就是看预测的p和实际y的匹配情况，这样我们可以得出匹配情况的百分比

可以看出89%

模型已经训练好了，然后我们就要将其可视化，第一要画出训练数据，第二我们要画出决策边界。

plotData(data, 'Exam 1 score', 'Exam 2 score', 'Pass', 'Failed')

def plotData(data, label_x, label_y, label_pos, label_neg, axes=None):

# 获得正负样本的下标(即哪些是正样本，哪些是负样本)

neg = data[:, 2] == 0

pos = data[:, 2] == 1

if axes == None:

axes = plt.gca()

axes.scatter(data[pos][:, 0], data[pos][:, 1], marker='+', c='k', s=60, linewidth=2, label=label_pos)

axes.scatter(data[neg][:, 0], data[neg][:, 1], c='y', s=60, label=label_neg)

axes.set_xlabel(label_x)

axes.set_ylabel(label_y)

axes.legend(frameon=True, fancybox=True)

这个就是画出所有的训练数据,neg为存储所有类别为0的，pos为存储所有类别为1的。

axes.scatter(data[pos][:, 0], data[pos][:, 1], marker='+', c='k', s=60, linewidth=2, label=label_pos)

这个是画出所有的正样本

axes.scatter(data[neg][:, 0], data[neg][:, 1], c='y', s=60, label=label_neg)

这个是画出所有的负样本

下面来画出决策的边界：

x1_min, x1_max = X[:,1].min(), X[:,1].max(),x2_min, x2_max = X[:,2].min(), X[:,2].max(),xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(res.x))h = h.reshape(xx1.shape)plt.contour(xx1, xx2, h, [0.5], linewidths=1, colors='b');plt.show()

h是对所有的网格数据进行预测，结果是0到1之间的值，然后我们通过contour来画数据，画h=0.5的那条线：

【0.01】时的效果是：

至此逻辑回归的全部代码

# -*- coding: utf-8 -*-import pandas as pdimport numpy as npimport matplotlib as mplimport matplotlib.pyplot as pltfrom scipy.optimize import minimizefrom sklearn.preprocessing import PolynomialFeaturesimport seaborn as snsdef loaddata(file, delimeter): data = np.loadtxt(file, delimiter=delimeter) print('Dimensions: ',data.shape) print(data[1:6,:]) return(data)def plotData(data, label_x, label_y, label_pos, label_neg, axes=None): # 获得正负样本的下标(即哪些是正样本，哪些是负样本) neg = data[:, 2] == 0 pos = data[:, 2] == 1 if axes == None: axes = plt.gca() axes.scatter(data[pos][:, 0], data[pos][:, 1], marker='+', c='k', s=60, linewidth=2, label=label_pos) axes.scatter(data[neg][:, 0], data[neg][:, 1], c='y', s=60, label=label_neg) axes.set_xlabel(label_x) axes.set_ylabel(label_y) axes.legend(frameon=True, fancybox=True)#定义sigmoid函数def sigmoid(z): return(1.0 / (1 + np.exp(-z)))def costFunction(theta, X, y): m = y.size h = sigmoid(X.dot(theta)) J = -1.0 * (1.0 / m) * (np.log(h).T.dot(y) + np.log(1 - h).T.dot(1 - y)) if np.isnan(J[0]): return (np.inf) return (J[0])def gradient(theta, X, y): m = y.size h = sigmoid(X.dot(theta.reshape(-1, 1))) grad = (1.0 / m) * X.T.dot(h - y) return (grad.flatten())def predict(theta, X, threshold=0.5): p = sigmoid(X.dot(theta.T)) >= threshold return(p.astype('int'))data = loaddata('data1.txt', ',')X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]y = np.c_[data[:,2]]initial_theta = np.zeros(X.shape[1])cost = costFunction(initial_theta, X, y)grad = gradient(initial_theta, X, y)print('Cost: \n', cost)print('Grad: \n', grad)res = minimize(costFunction, initial_theta, args=(X,y), method=None, jac=gradient, options={'maxiter':400})a=np.array([1, 10, 10])b=predict(res.x,a)print bp = predict(res.x, X)print('Train accuracy {}%'.format(100*sum(p == y.ravel())/p.size))plotData(data, 'Exam 1 score', 'Exam 2 score', 'Pass', 'Failed')x1_min, x1_max = X[:,1].min(), X[:,1].max(),x2_min, x2_max = X[:,2].min(), X[:,2].max(),xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(res.x))h = h.reshape(xx1.shape)plt.contour(xx1, xx2, h, [0.01], linewidths=1, colors='b');plt.show()

机器学习之逻辑回归向量化线性可分的