逻辑回归-Logistic Regression

1、前言

首先值得注意的是，逻辑回归（对数几率回归）虽然叫回归，但其实是个分类器。前面学习了线性回归，是用线性方程逼近一个高维平面，如果对y增加函数，则表示用线性方程去逼近一个高维曲面（大概可以这么理解吧），比如：

就是意图用线性方程逼近lny

2、问题引入

前面学习了线性回归，这里尝试用它来进行分类，引用周志华的《机器学习》，这部分写的很好：

不过我还是有个疑问，这里是阶跃函数不可导，因此要用sigmoid函数，但是在单层感知器中，用的是符号函数，应该也是不连续可导的，最后也能用梯度下降法求解（确实可以写出类似梯度的学习信号）。

几率一词并不陌生，这里就是尝试用线性方程去逼近事件发生的几率的对数。

平时我们可能会说明天下雨的几率很大，或者几率很小，其实就是对事件发生的可能性大小进行描述。生活中说概率还是比较小的，因为我们也不是明天下雨的概率是0.1还是0.8，一般都说几率大。

这里也指出，“逻辑回归”的翻译并不合适，一开始我也对这个名词感到疑惑

3、求解过程

其中有（周的似然项是加，按照他的似然项，我没有推出最后他得到的结果）：

这里已经把w和b合并一起为β（theta），求解β就是最大化损失函数L的过程。

将式（3.26）的p0消除并联立式（3.23）和式（3.24）代入（3.25）：（书上没有写该过程）

得到这样的导数，要求theta，直接用梯度下降法很容易求解，和线性回归差不多。值得注意的是，这里是求极大似然估计，是求最大值，为了保持下降的思维习惯，在梯度前面加个负号即可.不过必要性也不大，梯度下降是按照梯度的负方向行进，这里又有一个负号，因此最终还是theta+=梯度（）

4、代码实现

逻辑回归的代码与线性回归差不多，只是增加了sigmoid函数。为了验证效果，在二维平面上生成了两处随机点，并让他们的边缘部分重叠，如图：

sample_n = 1000
pos = np.random.normal(0, 0.1, sample_n)
pos = np.c_[pos, np.random.normal(0, 0.1, sample_n)]
y = np.ones([sample_n])

neg = np.random.normal(0, 0.1, sample_n) + 0.3
neg = np.c_[neg, np.random.normal(0, 0.1, sample_n) + 0.3]
y = np.r_[y, np.zeros(sample_n)]
y = y.reshape(y.shape[0], 1)
x = np.r_[pos, neg]
x = np.c_[x, np.ones(sample_n * 2)]  # 偏置项

这里用np.c_进行列合并，用np.r_进行行合并，最终生成样本为2*sample_n行，2列。

上面的推导中将b吸收入了w，所以这里x还要增加一列全是1的，而w将会有三个值

x = np.c_[x, np.ones(sample_n * 2)]  # 偏置项
w = np.zeros([x.shape[1], 1], dtype=np.float64)

下面使用批量梯度下降法来实现，至于SGD和BGD各自的优缺点在线性回归中已经分析并且全部做了代码实现，这里只用批量梯度下降法实现：

lr_init = 5
max_iter = 2000

for _iter in range(max_iter):
    lr = lr_init * (max_iter - _iter) / max_iter + 0.01
    y_ = x.dot(w)
    y_sig = sigmoid(y_)
    err = np.mean((y - y_sig) * x, axis=0)
    e2 = err.reshape(x.shape[1], 1)
    w += lr * e2
    if (_iter % 100 == 0) | (_iter == 1999):
        loss = np.mean(np.abs(y - sigmoid(x.dot(w))))
        print("iter", _iter, ":", w.reshape(w.shape[0]), loss)

为了加速收敛，学习率使用了线性下降的学习率，并设置最小学习率为0.01，后面的代码不难，不一一解释。注意的是err的shape和w并不一样，err的shape是一个但行向量，这里把它reshape和w一样的形状。最后计算一下损失，这里把损失定义为真实值与预测值的绝对值的均值

其中sigmoid函数为：

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

下面计算准确率：

# f1 = lambda a: 1 if a > 0.5 else 0
yo = [1 if b > 0.5 else 0 for b in sigmoid(x.dot(w))]
acc = 0.0
for j in range(x.shape[0]):
    if y[j] == yo[j]:
        acc += 1
acc /= x.shape[0]
print("acc rate:", acc)

其实这里可以不把预测值二值化，只是编程过程中显示方便，顺便做了一下这个工作。

最后是绘图：

# print(np.c_[sigmoid(x.dot(w)), yo])
plt.figure('aa')
plt.plot(pos[:, 0], pos[:, 1], 'bo')
plt.plot(neg[:, 0], neg[:, 1], 'yo')
k = -w[0] / w[1]  # 注意是0=w0*x0 + w1*x1 + w2,画图时候x1是纵轴
b = -w[2] / w[1]
_x = np.linspace(x[:, 0].min(), x[:, 0].max(), 50)
plt.plot(_x, _x*k + b, "r")
plt.show()

画直线时候用了样本x0的最小值和最大值，来让显示更加美观合理。线性方程斜率和偏置的求解不再赘述。

完整代码如下：

def logistic_regression():  # x_data, w_ideal, bias, noise_sigma
    sample_n = 1000
    pos = np.random.normal(0, 0.1, sample_n)
    pos = np.c_[pos, np.random.normal(0, 0.1, sample_n)]
    y = np.ones([sample_n])

    neg = np.random.normal(0, 0.1, sample_n) + 0.3
    neg = np.c_[neg, np.random.normal(0, 0.1, sample_n) + 0.3]
    y = np.r_[y, np.zeros(sample_n)]
    y = y.reshape(y.shape[0], 1)
    x = np.r_[pos, neg]

    x = np.c_[x, np.ones(sample_n * 2)]  # 偏置项
    w = np.zeros([x.shape[1], 1], dtype=np.float64)

    lr_init = 5
    max_iter = 2000
    # print(x)
    # print(y.shape)
    # print(w)

    for _iter in range(max_iter):
        lr = lr_init * (max_iter - _iter) / max_iter + 0.01
        y_ = x.dot(w)
        y_sig = sigmoid(y_)
        err = np.mean((y - y_sig) * x, axis=0)
        e2 = err.reshape(x.shape[1], 1)
        w += lr * e2
        if (_iter % 100 == 0) | (_iter == 1999):
            loss = np.mean(np.abs(y - sigmoid(x.dot(w))))
            print("iter", _iter, ":", w.reshape(w.shape[0]), loss)

    # f1 = lambda a: 1 if a > 0.5 else 0
    yo = [1 if b > 0.5 else 0 for b in sigmoid(x.dot(w))]
    acc = 0.0
    for j in range(x.shape[0]):
        if y[j] == yo[j]:
            acc += 1
    acc /= x.shape[0]
    print("acc rate:", acc)

    # print(np.c_[sigmoid(x.dot(w)), yo])
    plt.figure('aa')
    plt.plot(pos[:, 0], pos[:, 1], 'bo')
    plt.plot(neg[:, 0], neg[:, 1], 'yo')
    k = -w[0] / w[1]  # 注意是0=w0*x0 + w1*x1 + w2,画图时候x1是纵轴
    b = -w[2] / w[1]
    _x = np.linspace(x[:, 0].min(), x[:, 0].max(), 50)
    plt.plot(_x, _x*k + b, "r")
    plt.show()

运行效果：

iter 0 : [-0.37836782 -0.37063006  0.        ] 0.472136956147
iter 100 : [-9.68396276 -9.31203219  2.79515917] 0.0972483530724
iter 200 : [-12.3732208  -11.77084137   3.57636014] 0.0696249130477
iter 300 : [-14.05242298 -13.25502332   4.05825691] 0.0585232694782
iter 400 : [-15.27801653 -14.30871006   4.40623915] 0.0523445352026
iter 500 : [-16.23755579 -15.11479075   4.6761783 ] 0.048357392951
iter 600 : [-17.01856995 -15.75819429   4.89415953] 0.0455577543677
iter 700 : [-17.66944873 -16.28555085   5.07458808] 0.0434832090325
iter 800 : [-18.21992199 -16.72527149   5.22629481] 0.041888970143
iter 900 : [-18.68961704 -17.09593876   5.3550923 ] 0.0406325122006
iter 1000 : [-19.09217625 -17.41034163   5.4650066 ] 0.0396249222964
iter 1100 : [-19.43744718 -17.67761876   5.55893285] 0.0388077837751
iter 1200 : [-19.73273803 -17.90448472   5.63901127] 0.0381411686555
iter 1300 : [-19.98358021 -18.09597288   5.70685554] 0.037596954919
iter 1400 : [-20.19421293 -18.25590598   5.76369782] 0.0371548985894
iter 1500 : [-20.3679022  -18.38720493   5.81048421] 0.0368002225101
iter 1600 : [-20.5071564  -18.49209662   5.84793912] 0.0365220852894
iter 1700 : [-20.61387433 -18.57225617   5.87660948] 0.0363125854021
iter 1800 : [-20.68944732 -18.6289047    5.89689518] 0.0361661053985
iter 1900 : [-20.73482879 -18.66287567   5.90906982] 0.0360788825377
iter 1999 : [-20.7505675  -18.67464901   5.91329088] 0.036048761143
acc rate: 0.985