逻辑回归实现

逻辑斯谛分布

• 定义(逻辑斯谛分布):设X是连续随机变量,X服从逻辑斯谛分布是指X具有下列分布函数和密度函数
  $$F(x)=P(X\le x)=\frac{1}{1+e^{-(x-\mu )/\gamma }}$$

$$f(x)=F^{{}^{\prime }}(x)=\frac{e^{-(x-\mu )/\gamma }}{\gamma {(1+e^{-(x-\mu )/\gamma })}^2}$$

其中$\mu$为位置参数,$\gamma$为形状参数

很多人对这个分布没有直观的理解,话不多说 我们直接上图 这个分布如图所示

import matplotlib.pyplot as plt
import numpy as np
mu=0
gamma=[0.5,1,2]
x_range=np.arange(-10,10,0.1)
fig,(ax1,ax2)=plt.subplots(1,2)
for g in gamma:
    ax1.plot(x_range,1/(1+np.exp(-(x_range-mu)/g)),label=f"$\mu$={mu},$\gamma$={g}")
    ax2.plot(x_range,(np.exp(-(x_range-mu)/g))/(g*(1+np.exp(-(x_range-mu)/g))**2),label=f"$\mu$={mu},$\gamma$={g}")
ax1.legend()
ax2.legend()
ax1.set_title("F(x)")
ax2.set_title("f(x)")
ax1.grid()
ax2.grid()
plt.subplots_adjust(left=0,right=2,top=1,hspace=0.01)

二项逻辑斯谛回归模型

二项逻辑斯谛回归模型是一种分类模型,由条件概率分布P(Y|X)表示,形式为参数化的逻辑斯谛分布,随机变量X取值实数,随机变量Y取值为0或为1,通过监督学习的方式估计模型参数

• 定义:(逻辑斯谛回归模型)

$$P(Y=1|x)=\frac{exp(w\cdot x+b)}{1+exp(w\cdot x+b)}$$

$$P(Y=0|x)=\frac{1}{1+exp(w\cdot x+b)}$$

$x\in R^n$是输入,Y ∈ { 0 , 1 } \in\{0,1\} ∈{0,1}是输出,$w\in R^n$和b$\in R$是参数,$w$称为权重,b称为偏置,$w\cdot x$为$w$和$x$的内积
为了考察逻辑斯谛回归模型的特点,一个事件的几率(odds)是指改时间发生的概率与该事件不发生的概率的比值,如果事件发生的概率是p,那么该事件的几率是$\frac{p}{1-p}$,该事件的对数几率(log odds)或logit函数是
$$logit(p)=log\frac{p}{1-p}$$
对于逻辑斯谛回归而言 由(3)和(4)可得
$$\begin{gathered} log\frac{P(Y=1|x)}{1-P(Y=1|x)} \\ =log\frac{P(Y=1|x)}{P(Y=0|x)}=logexp(w\cdot x)=w\cdot x \end{gathered}$$

在逻辑斯谛回归模型当中,输出Y=1的对数几率是输入x的线性函数或者说输出Y=1的对数几率是由输入线性函数表示的模型,即逻辑斯谛模型回归模型,换一个角度来看,考虑对输入x进行分类的线性函数$w\cdot x$其值域为实数域 ,注意这里的x$\in R^{n+1}w\in R^{n+1}$,通过逻辑斯谛回归模型定义式将线性函数$w\cdot x$转化为概率
$$P(Y=1|x)=\frac{exp(w\cdot x)}{1+exp(w\cdot x)}$$

若线性函数的值越接近正无穷概率值就为1,接近于负无穷概率值为0,这样的模型我们成为逻辑斯谛回归模型

模型参数估计

逻辑斯谛回归模型学习时,对于给定的训练数据 T = { ( x 1 , y 1 ) , ( x 2 , y 2 ) . . . ( x n , y n ) } T=\{(x_1,y_1),(x_2,y_2)...(x_n,y_n)\} T={(x1​,y1​),(x2​,y2​)...(xn​,yn​)},其中$x_i\in R^n,y_i\in 0,1\}$可以应用极大似然估计法估计模型参数

设:
$$P(Y=1|x)=\pi (x),P(Y=0|x)=1-\pi (x)$$

似然函数为:

$$\prod _{i=1}^n{[\pi (x_i)]}^{y_i}{[1-\pi (x_i)]}^{1-y_i}$$

对数似然函数为:
$$L(w)=\sum _{i=1}^n[y_{i}log\pi (x_i)+(1-y_i)log(1-\pi (x_i))]$$

对L(w)求极大值,得到w的估计值(通常用梯度下降法和拟牛顿法得到最佳参数$\hat{w}$)
到了这里我们发现 对数似然函数就等于损失函数的相反数
因为我们发现如果$y_i=1$ ,只有$y_{i}log(\pi x(i))$这一项生效,结果为0,若$y_i=0$则只有$(1-y_i)log(1-\pi (x_i)这一项生效,结果仍然为0$,另外我们用
sigmoid函数代表所谓的逻辑斯谛分布
$$Sigmoid(x)=\frac{1}{1+e^{-x}}$$

损失函数求解

定义函数$J(w)=\frac{-L(w)}{n}$这里我们取平均损失
随机初始化 w = { w 1 , w 2 , . . . w n } 和 学 习 率 η w=\{w_1,w_2,...w_n\}和学习率\eta w={w1​,w2​,...wn​}和学习率η
$$\frac{\partial J}{\partial w_j}=\frac{\partial J}{\partial L}\times \frac{\partial L}{\partial w_j}$$
$$=\frac{\partial J}{\partial L}\times \frac{\partial L}{\partial \pi (x_i)}\times \frac{\partial \pi (xi)}{\partial (w\cdot x)}\times \frac{\partial (w\cdot x)}{\partial w_j}$$
$$=\frac{1}{n}\sum _{i=1}^n(\frac{-y_i}{\pi (x_i)}+\frac{1-y_i}{1-\pi (x_i)})(\pi (x_i)(1-\pi (x_i)))x_j$$
$$=\frac{1}{n}\sum _{i=1}^n(\frac{-y_i}{\pi (x_i)}+\frac{1-y_i}{1-\pi (x_i)})(\pi (x_i)(1-\pi (x_i)))x_j$$
$$=\frac{1}{n}\sum _{i=1}^n(\pi (x_i)-y_i)x_i^j$$

其中$\pi (x_i)=sigmoid(w\cdot x),x_i^j代表第i个样本第j个特征$

代码实现

数据加载

数据集链接 提取码p7gq

import numpy as np
import matplotlib.pyplot as plt
#加载数据集
with open("testData.txt",'r') as reader:
    txt=reader.read()
data=np.array([list(map(eval,item.split())) for item in 
txt.split("\n")])
labels=data[:,2].copy()#获取标签y
y1_index=data[:,2]==0
y2_index=data[:,2]==1
data[:,2]=1#此处其实是把b当成特征常量1
plt.scatter(data[:,0][y1_index],data[:,1][y1_index])
plt.scatter(data[:,0][y2_index],data[:,1][y2_index])

数据分布如图: 我们可以发现这其实不是一个比较好的线性分布

训练过程

def g(x):#sigmoid 函数
    return 1/(1+np.exp(-x))
w=np.array([7.0,6.0,5.0])
epoches=500#迭代500轮
lr=1#学习率设置为1
nums=len(labels)#获取样本总数
cost_records=[]#记录损失值
for epoche in range(epoches):
    dldw1,dldw2,dldw3=0,0,0
    for num in range(nums):
        dldw1+=(g(w@data[num])-labels[num])*data[num][0]
        dldw2+=(g(w@data[num])-labels[num])*data[num][1]
        dldw3+=(g(w@data[num])-labels[num])*data[num][2]
    cost=np.sum(np.abs(g(data@w.T)-labels))/nums
    cost_records.append(cost)
    lr_update=lr*0.5**(epoche/50)#这里采用指数学习率衰减
    w[0]=w[0]-lr_update*dldw1/nums
    w[1]=w[1]-lr_update*dldw2/nums
    w[2]=w[2]-lr_update*dldw3/nums

验证算法

$$\begin{gathered} \hat{w}\cdot x=0 \\ =w_1{x}_1+w_2{x}_2+w_3(b) \end{gathered}$$
我们随机生成坐标来拟合回归直线

testx=np.arange(-4,4,0.5)
testy=(-w[2]-w[0]*testx)/w[1]
plt.scatter(data[:,0][y1_index],data[:,1][y1_index])
plt.scatter(data[:,0][y2_index],data[:,1][y2_index])
plt.plot(testx,testy)

效果还是很不错的,其中有几个误分类点,这是因为原始数据的分布不是完美的线性

Loss损失:

plt.plot(list(range(epoches)),cost_records)
plt.title("$Learning Rate: \eta $")

如有错误,敬请指正