逻辑回归手动实现（logistic regression）

受到生物神经网络结构的启发，模拟生物神经网络结构设计了以下的算法

假设X是一个2维的特征
h ( X ; w ) = s i g m o i d ( x 2 ∗ w 2 + x 1 ∗ w 1 + w 0 ) = s i g m o i d ( X w ) h(X;w) = sigmoid(x_2*w_2+x_1*w_1+w_0)=sigmoid(Xw) h(X;w)=sigmoid(x2​∗w2​+x1​∗w1​+w0​)=sigmoid(Xw)
其中
     σ ( z ) = s i g m o i d ( z ) = 1 1 + e − z σ(z) = sigmoid(z)=\frac{1}{1+e^{-z}} σ(z)=sigmoid(z)=1+e−z1​

sigmoid函数图像如下

输出范围[0,1],可以把输出理解为Y=1的概率概率>0.5 输出标签1，否则输出0

已知输入特征Xi ， 预测输出标签Yi=1 的概率为：
P ( Y i = 1 ∣ X i ; w ) = h ( X i ; w ) = σ ( X i w ) P( Y_i=1|X_i;w) = h(X_i;w) = σ(X_iw) P(Yi​=1∣Xi​;w)=h(Xi​;w)=σ(Xi​w)

已知输入特征Xi ， 预测输出标签Yi=0 的概率为：
P ( Y i = 0 ∣ X i ; w ) = 1 − h ( X i ; w ) P( Y_i=0|X_i;w) = 1-h(X_i;w) P(Yi​=0∣Xi​;w)=1−h(Xi​;w)

已知输入特征Xi以及对应的标签Yi
则预测准确的概率为：

P ( Y i , X i ; w ) = [ h ( X i ; w ) ] Y i [ 1 − h ( X i ; w ) ] 1 − Y i P( Y_i,X_i;w) = [h(X_i;w)]^{Y_i}[1-h(X_i;w)]^{1-Y_i} P(Yi​,Xi​;w)=[h(Xi​;w)]Yi​[1−h(Xi​;w)]1−Yi​
预测的越准 上述函数就越大
举例说明：
假设真实标签为1，预测概率为0.8；则上述式子为 0. 8 1 ∗ ( 1 − 0.8 ) 0 = 0.8 0.8^1*(1-0.8)^0=0.8 0.81∗(1−0.8)0=0.8
假设真实标签为0，预测概率为0.8；则上述式子为 0. 8 0 ∗ ( 1 − 0.8 ) 1 = 0.2 0.8^0*(1-0.8)^1=0.2 0.80∗(1−0.8)1=0.2
可以看出上述函数就越大，预测的越准；

根据已知数据，寻找最优参数 w，就是令

一般对连积形式，我们都会取log形式
取log写成似然函数的形式

这个公式看上去是否有那么一丢丢熟悉，前面我们在朴素贝叶斯分类里也出现了类似的公式 P ∗ l o g ( P ) P*log(P) P∗log(P) 还记得这个公式嘛，这个公式就是信息熵的公式

信息熵的值越大，说明系统越混乱；信息熵越小系统越确定,分类也就越准确

但这里的公式与前面的信息熵公式看起来又有那么一点点不同，前面是自己和自己的log相乘，这里多了一项，这个公式叫做交叉熵，作用和信息熵差不多

将上面的公式再加一个“-”号上式就变成了了熵的形式（交叉熵）

熵越小系统越确定,分类也就越准确

逻辑回归的最终目标变为：
已知训练数据X以及相应标签Y求解参数w使 交叉熵L(w) 最小；

不同于线性回归，这里的公式不能直接求出闭式解，只能通过别的方式来求出最优解（尽可能的最优）

这里就需要引出梯度下降的概念

简单来说就是对原函数求它的一阶导数，沿着当前w的切线的反方向移动一小步，然后再次求移动后的w处的切线，沿着它的反方向移动一小步，重复上述过程，直到达到最优点。

同理，当w出现在右侧时

每次都沿着梯度的反方向前进一小段，这个一小段是由公式中的 α 确定的，α 的值越小，每次前进的步伐就小一点， α 的值越大，每次前进的就多一点；
所以 α 的值很关键，如果 α 的值过大，就会出现反复横跳的现象

反复横跳，可能会导致函数最终不会收敛，导致不能得到准确的预测，所以 α 的值一般都设置的比较小，0.001-0.01之间；具体可以通过调试得到；

下面给出公式的化简过程，以及转换为矩阵形式，感兴趣的可以看一下，不感兴趣其实也没多大影响，现在很多包都可以直接调用。

引入相关包，并定义sigmoid函数

import numpy as np
import matplotlib.pyplot as plt

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

权重更新函数

# datas NxD
# labs Nx1
# w    Dx1

def weight_update(datas,labs,w,alpha=0.01):
    z = np.dot(datas,w) # Nx1
    h = sigmoid(z)        # Nx1
    Error = labs-h        # Nx1 
    w = w + alpha*np.dot(datas.T,Error)
    return w

训练函数

def train_LR(datas,labs,n_epoch=2,alpha=0.005):
    N,D = np.shape(datas)   
    w = np.ones([D,1])  # Dx1
    # 进行n_epoch轮迭代
    for i in range(n_epoch):
        w = weight_update(datas,labs,w,alpha)
        error_rate=test_accuracy(datas,labs,w)
        print("epoch %d error %.3f%%"%(i,error_rate*100))
    return w

计算准确率

def test_accuracy(datas,labs,w):
    N,D = np.shape(datas)
    z = np.dot(datas,w) # Nx1
    h = sigmoid(z)        # Nx1
    lab_det = (h>0.5).astype(np.float)
    error_rate=np.sum(np.abs(labs-lab_det))/N
    return error_rate

画出决策边界

def draw_desion_line(datas,labs,w,name="0.jpg"):
    dic_colors={0:(.8,0,0),1:(0,.8,0)}
  
    # 画数据点
    for i in range(2):
        index = np.where(labs==i)[0]
        sub_datas = datas[index]
        plt.scatter(sub_datas[:,1],sub_datas[:,2],s=16.,color=dic_colors[i])
    
    # 画判决线
    min_x = np.min(datas[:,1])
    max_x = np.max(datas[:,1])
    w = w[:,0]
    x = np.arange(min_x,max_x,0.01)
    y = -(x*w[1]+w[0])/w[2]
    plt.plot(x,y)
    
    plt.savefig(name)

加载数据集

def load_dataset(file):    
    with open(file,"r",encoding="utf-8") as f:
        lines = f.read().splitlines()
    
    # 取 lab 维度为 N x 1
    labs = [line.split("\t")[-1] for line in lines]
    labs = np.array(labs).astype(np.float32)
    labs= np.expand_dims(labs,axis=-1) # Nx1
    
    # 取数据 增加 一维全是1的特征
    datas = [line.split("\t")[:-1] for line in lines]
    datas = np.array(datas).astype(np.float32)
    N,D = np.shape(datas)
    # 增加一个维度
    datas = np.c_[np.ones([N,1]),datas]
    return datas,labs

主函数

if __name__ == "__main__":
    ''' 实验1 基础测试数据'''
    # 加载数据
    file = "testset.txt"
    datas,labs = load_dataset(file)
    
    weights = train_LR(datas,labs,alpha=0.001,n_epoch=800)
    print(weights)
    draw_desion_line(datas,labs,weights,name="test_1.jpg")
    

运行结果

上面程序中使用的训练函数定义如下

def train_LR(datas,labs,n_epoch=2,alpha=0.005):
    N,D = np.shape(datas)   
    w = np.ones([D,1])  # Dx1
    # 进行n_epoch轮迭代
    for i in range(n_epoch):
        w = weight_update(datas,labs,w,alpha)
        error_rate=test_accuracy(datas,labs,w)
        print("epoch %d error %.3f%%"%(i,error_rate*100))
    return w

每一次都对数据集里面所有的数据求梯度，这样当数据集很大时，计算的效率就十分的缓慢，还可能会出现过拟合的现象；

为了解决这个问题，就可以使用随机梯度下降的方法，每次只从数据集里面取出一个bath，对bath里面的数据求梯度来更新权重w，实验表明，随机梯度下降的算法不仅可以提高运算效率，往往还能获得更精确的预测结果；

将数据分成小的Batch，每个Batch更新一次，可以增加训练的速度和稳定性

# 随机梯度下降
def train_LR_batch(datas,labs,batchsize,n_epoch=2,alpha=0.005):
    N,D = np.shape(datas)
    # weight 初始化
    w = np.ones([D,1])  # Dx1
    N_batch = N//batchsize
    for i in range(n_epoch):
        # 数据打乱
        rand_index  = np.random.permutation(N).tolist()
        # 每个batch 更新一下weight
        for j in range(N_batch):
            # alpha = 4.0/(i+j+1) +0.01
            index = rand_index[j*batchsize:(j+1)*batchsize]
            batch_datas = datas[index]
            batch_labs = labs[index]
            w=weight_update(batch_datas,batch_labs,w,alpha)
    
        error = test_accuracy(datas,labs,w)
        print("epoch %d  error  %.2f%%"%(i,error*100))
    return w

将上面代码主函数中的训练函数改为随机梯度下降即可；

testset.txt

-0.017612	14.053064	0
-1.395634	4.662541	1
-0.752157	6.538620	0
-1.322371	7.152853	0
0.423363	11.054677	0
0.406704	7.067335	1
0.667394	12.741452	0
-2.460150	6.866805	1
0.569411	9.548755	0
-0.026632	10.427743	0
0.850433	6.920334	1
1.347183	13.175500	0
1.176813	3.167020	1
-1.781871	9.097953	0
-0.566606	5.749003	1
0.931635	1.589505	1
-0.024205	6.151823	1
-0.036453	2.690988	1
-0.196949	0.444165	1
1.014459	5.754399	1
1.985298	3.230619	1
-1.693453	-0.557540	1
-0.576525	11.778922	0
-0.346811	-1.678730	1
-2.124484	2.672471	1
1.217916	9.597015	0
-0.733928	9.098687	0
-3.642001	-1.618087	1
0.315985	3.523953	1
1.416614	9.619232	0
-0.386323	3.989286	1
0.556921	8.294984	1
1.224863	11.587360	0
-1.347803	-2.406051	1
1.196604	4.951851	1
0.275221	9.543647	0
0.470575	9.332488	0
-1.889567	9.542662	0
-1.527893	12.150579	0
-1.185247	11.309318	0
-0.445678	3.297303	1
1.042222	6.105155	1
-0.618787	10.320986	0
1.152083	0.548467	1
0.828534	2.676045	1
-1.237728	10.549033	0
-0.683565	-2.166125	1
0.229456	5.921938	1
-0.959885	11.555336	0
0.492911	10.993324	0
0.184992	8.721488	0
-0.355715	10.325976	0
-0.397822	8.058397	0
0.824839	13.730343	0
1.507278	5.027866	1
0.099671	6.835839	1
-0.344008	10.717485	0
1.785928	7.718645	1
-0.918801	11.560217	0
-0.364009	4.747300	1
-0.841722	4.119083	1
0.490426	1.960539	1
-0.007194	9.075792	0
0.356107	12.447863	0
0.342578	12.281162	0
-0.810823	-1.466018	1
2.530777	6.476801	1
1.296683	11.607559	0
0.475487	12.040035	0
-0.783277	11.009725	0
0.074798	11.023650	0
-1.337472	0.468339	1
-0.102781	13.763651	0
-0.147324	2.874846	1
0.518389	9.887035	0
1.015399	7.571882	0
-1.658086	-0.027255	1
1.319944	2.171228	1
2.056216	5.019981	1
-0.851633	4.375691	1
-1.510047	6.061992	0
-1.076637	-3.181888	1
1.821096	10.283990	0
3.010150	8.401766	1
-1.099458	1.688274	1
-0.834872	-1.733869	1
-0.846637	3.849075	1
1.400102	12.628781	0
1.752842	5.468166	1
0.078557	0.059736	1
0.089392	-0.715300	1
1.825662	12.693808	0
0.197445	9.744638	0
0.126117	0.922311	1
-0.679797	1.220530	1
0.677983	2.556666	1
0.761349	10.693862	0
-2.168791	0.143632	1
1.388610	9.341997	0
0.317029	14.739025	0