机器学习之逻辑回归

上一章学习了一般线性回归，现在我们来看一下逻辑回归到底是什么鬼？

其实从这点上我认为，逻辑回归已经不属于什么回归了，而直接属于一种分类问题，但分类，为什么还叫回归呢，因为我们的逻辑回归使用的函数是这样的：

如果预测大于0.5，就判断为正，否则为负，对于回归，我们必不可少的函数就是损失函数，一般线性回归的函数如下：

而逻辑回归的损失函数是这样的：

这不是一个凸函数，因此很难找到一个最优解，这时候，数学家发挥作用了，他们用他们的智慧，研究出了下面的模型：

也就是这样的：

让我们来实现一下逻辑回归到底是什么样的

首先读入数据：

import matplotlib.pyplot as plt
from scipy.optimize import minimize

def 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);

那我们来读取数据吧：

data = loaddata(‘data1.txt‘,‘,‘)

我们来看一下数据是什么样的？

X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]
y = np.c_[data[:,2]]
plotData(data, ‘Exam 1 score‘, ‘Exam 2 score‘, ‘Pass‘, ‘Fail‘)

因此我们的任务就是找到一个分界点很好的将这两类数据分开。

定义逻辑回归

def sigmoid(z):
    return(1 / (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()

我们先来看看初始值的误差和梯度：

initial_theta = np.zeros(X.shape[1])
cost = costFunction(initial_theta,X,y) #0.69314718055994529
grad = gradient(initial_theta,X,y)     #[ -0.1 , -12.00921659, -11.26284221]

这样我们就来最小化这个函数：

res = minimize(costFunction, initial_theta, args=(X,y), jac=gradient, options={‘maxiter‘:400})  #这里使用了scipy中的库，不懂的同学可以看看官方文档，很简单的。

那我们的预测函数时这样的：

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

我们将前面的数据集画出边界：

plt.scatter(45, 85, s=60, c=‘r‘, marker=‘v‘, label=‘(45, 85)‘)
plotData(data, ‘Exam 1 score‘, ‘Exam 2 score‘, ‘Admitted‘, ‘Not admitted‘)
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‘);

带有正则化的逻辑回归

前面我们学习了一般的逻辑回归，没有带有正则化参数，我们来看一下正则化怎么弄？

data2 = loaddata(‘data2.txt‘,‘,‘)
X = data2[:,0:2]
y = np.c_[data2[:,2]]

plotData(data2, ‘Microchip Test 1‘, ‘Microchip Test 2‘, ‘y = 1‘, ‘y = 0‘)
plt.show()

我们需要找一个边界很好地分类这些点。可以看见，点与点之间并不是严格分开的，我们自己定义一个高阶函数来拟合这个数据集。

from sklearn.preprocessing import PolynomialFeatures

poly = PolynomialFeatures(6)
XX = poly.fit_transform(data2[:,0:2])

def costFunctionReg(theta, reg, *args):
    m = y.size
    h = sigmoid(XX.dot(theta))

    J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y)) + (reg/(2.0*m))*np.sum(np.square(theta[1:]))

    if np.isnan(J[0]):
        return(np.inf)
    return(J[0])

def gradientReg(theta, reg, *args):
    m = y.size
    h = sigmoid(XX.dot(theta.reshape(-1,1)))

    grad = (1.0/m)*XX.T.dot(h-y) + (reg/m)*np.r_[[[0]],theta[1:].reshape(-1,1)]

    return(grad.flatten())

那我们来看看不同的正则化系数对结果的影响：

initial_theta = np.zeros(XX.shape[1])
costFunctionReg(initial_theta, 1, XX, y)

fig, axes = plt.subplots(1,3, sharey = True, figsize=(17,5))

# 决策边界，咱们分别来看看正则化系数lambda太大太小分别会出现什么情况
# Lambda = 0 : 就是没有正则化，这样的话，就过拟合咯
# Lambda = 1 : 这才是正确的打开方式
# Lambda = 100 : 正则化项太激进，导致基本就没拟合出决策边界

for i, C in enumerate([0.0, 1.0, 100.0]):
    # 最优化 costFunctionReg
    res2 = minimize(costFunctionReg, initial_theta, args=(C, XX, y), jac=gradientReg, options={‘maxiter‘:3000})

    # 准确率
    accuracy = 100.0*sum(predict(res2.x, XX) == y.ravel())/y.size    

    # 对X,y的散列绘图
    plotData(data2, ‘Microchip Test 1‘, ‘Microchip Test 2‘, ‘y = 1‘, ‘y = 0‘, axes.flatten()[i])

    # 画出决策边界
    x1_min, x1_max = X[:,0].min(), X[:,0].max(),
    x2_min, x2_max = X[:,1].min(), X[:,1].max(),
    xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))
    h = sigmoid(poly.fit_transform(np.c_[xx1.ravel(), xx2.ravel()]).dot(res2.x))
    h = h.reshape(xx1.shape)
    axes.flatten()[i].contour(xx1, xx2, h, [0.5], linewidths=1, colors=‘g‘);       
    axes.flatten()[i].set_title(‘Train accuracy {}% with Lambda = {}‘.format(np.round(accuracy, decimals=2), C))

到这里，回归已经学习完毕，在下还在学习阶段，如有不对之处，请吐槽！