【机器学习】PyTorch手动实现Logistic回归算法

本算法算是本人在学习机器学习过程中的入门算法之一。

该算法实现了纯手工打造，每一行代码都包含原创作者的心血，是算法中的劳斯莱斯，向该算法的原创作者致敬！

参考地址：点击打开

计算较为繁琐，需要用到sigmoid函数和梯度下降算法，算法的步骤主要如下：

1. 二项分布概率公式
2. 最大似然估计法
3. 求导数
4. 梯度下降算法和优化

上图的主要过程是对交叉熵里的w求导，手写的拍照，可能有些杂乱，而且略有错误，下面再次手动推导一下。

交叉熵是概率统计里一种描述预测差异的方法，主要思路是单分类任务的结果服从二项分布，而且其出现一次的概率为：

$$P(x=k)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)P_{}^k(1-P{)}_{}^{n-k}$$

其总的概率分布为：
$$\sum _{i=1}^{n}P(x=i)=\sum _{i=i}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)P_{}^i(1-P{)}_{}^{n-i}=1$$

因为二分类任务的输出一般只有两种（0和1），用y 和1-y来表示，因此产生的交叉熵为：$-log{p}_i^y(1-p_i{)}^{1-y}$，对此式取对数并化简：
$$\begin{gathered} -log(p_i^y(1-p_i{)}^{1-y})=-[ylog{p}_i+(1-y)log(1-p_i)]= \\ -[ylog{p}_i+log(1-p_i)-ylog(1-p_i)]= \\ -[ylog\frac{p_i}{1-p_i}+log(1-p_i)] \end{gathered}$$

对w求偏导：
$$\begin{gathered} \frac{\partial f}{\partial w}=-\frac{\partial log(p_i^y(1-p_i{)}^{1-y})}{\partial w}= \\ -\frac{\partial [ylog\frac{p_i}{1-p_i}+log(1-p_i)]}{\partial w}= \\ -\frac{\partial [ylog\frac{\frac{1}{1+e^{-wx}}}{1-\frac{1}{1+e^{-wx}}}+log(1-\frac{1}{1+e^{-wx}})]}{\partial w}= \\ -\frac{\partial [ylog\frac{\frac{1}{1+e^{-wx}}}{\frac{e^{-wx}}{1+e^{-wx}}}+log(\frac{e^{-wx}}{1+e^{-wx}})]}{\partial w}= \\ -\frac{\partial [ylog\frac{1}{e^{-wx}}+log(\frac{e^{-wx}}{1+e^{-wx}})]}{\partial w}= \\ -\frac{\partial [ywx+log(\frac{e^{-wx}}{1+e^{-wx}})]}{\partial w}= \\ -yx-\frac{\partial [log(\frac{e^{-wx}}{1+e^{-wx}})]}{\partial w}= \\ -yx-\frac{1+e^{-wx}}{e^{-wx}}\frac{-e^{-wx}(-x)e^{-wx}+(1+e^{-wx})e^{-wx}(-x)}{[1+e^{-wx}{]}^2}= \\ -yx-\frac{1+e^{-wx}}{e^{-wx}}\frac{e^{-2wx}(x)-e^{-wx}x-x{e}^{-2wx}}{[1+e^{-wx}{]}^2}= \\ -yx-\frac{-x}{1+e^{-wx}}=(\hat{y}-y)x \end{gathered}$$

上图中，y和$\hat{y}$的顺序能不能交换呢？测试是不可以的，为什么呢？此处还有存疑，有时间再好好想想。

进一步的阐述说明：

1. sigmoid函数没有偏执bias参数。
2. 上面导数推导图中的连乘号不准确，实际是表达了预测结果与实际结果两个矩阵向量的点积运算。
3. 算法测度损失使用交叉熵。从熵的性质可以知道，交叉熵是熵的最大值，使用交叉熵的主要原因是便于计算吗？
4. 二分类问题中，y 和 $\hat{y}$的概率和是1，次数之和是n，而不是y 的值为1，$\hat{y}$的值为0；如果是这样的话，就没法求导了。
5. 划线模块draw_desion_line的原理是怎样的呢？在模块load_dataset中对x增加一个维度的原因是什么呢？此问题是二分类问题画图的核心问题之一 ，初学者一般都会提出这样一个问题，这个问题的实质请参考如下网址的解释说明：https://blog.csdn.net/xq151750111/article/details/121871627，以及：https://scipython.com/blog/plotting-the-decision-boundary-of-a-logistic-regression-model/；https://stats.stackexchange.com/questions/202668/calculation-of-decision-boundaries-with-perceptron。回到第二个问题，测试集也可以将第一列设置为0，不影响最后的画图。

代码：

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_batch(datas,labs,batchsize=80,n_epoch=2,alpha=0.005):
    
    print("epoch:%d,alpha:%.8f batchsize:%d"%(n_epoch,alpha,batchsize))
    
    N,D = np.shape(datas)
    
    print("N:%f D:%f"% (N,D))
    
    # weight 初始化
    w = np.ones([D,1])  # Dx1
    N_batch = N//batchsize  #取整数
    
    #print("n:%d d:%d batchsize:%d"%(N,D,batchsize))
    
    for i in range(n_epoch):
        
        # 数据打乱
        rand_index  = np.random.permutation(N).tolist()
        # 每个batch 更新一下weight
        for j in range(N_batch):
            #print("i:%d j:%d N_batch:%d\r\n"%(i,j,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


def train_LR(datas,labs,n_epoch=2,alpha=0.005):
    N,D = np.shape(datas)   
    
    print("N:%f D:%f"% (N,D))
    
    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)}
  
    #print(labs)
    
    # 画数据点
    for i in range(2):
        #print(i)
        index = np.where(labs==i)[0]

        sub_datas = datas[index]
        #print(sub_datas)
        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)
    
    #print(x)
    
    y = -(x*w[1]+w[0])/w[2]
    
    #y = -(x*w[2]+w[1])/w[0]
    
    plt.plot(x,y)
    
    #print(w)
    
    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
    
    #rint(labs)

    # 取数据 增加 一维全是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]
    
    #print(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)
    
    weights = train_LR_batch(datas,labs,batchsize=80,alpha=0.001,n_epoch=800)
    print(weights)
    draw_desion_line(datas,labs,weights,name="test_1.jpg")
    

训练需要的数据集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

anaconda jupyter notebook中的运行结果截图：