logistics regression逻辑回归

问题引入

已知下面表格中的数据符合某种线性关系，对于预测某个人属于哪一个类别（1代表成人，0代表未成年人），存在下面的表格，现在需要利用一个人的视力身高以及体重信息，对这个人进行分类。

视力${\text{x}}_1$	身高${\text{x}}_2$	体重${\text{x}}_3$	类别$\text{y}$
3.8	180	85	1
3.7	186	81	1
3.8	190	100	1
5.2	160	50	0
5.1	166	70	0
5.0	170	60	0

---

符号说明

${\text{x}}_i=(1,x_{i1},x_{i2},...,x_{in}{)}^T$，是列向量，其中$x_{i1}\cdots x_{i_n}$表示表格中的非类别行。
${\text{x}}_j=(x_{1j},x_{2j},...,x_{mj}{)}^T$，是列向量，表示表格中的非类别列。
$\text{x}$，是个mxn的矩阵，表示表格中不包括类别的矩阵。
$\text{y}=(y_1,y_2,...,y_m{)}^T$，是列向量，表示表格中的类别所在列。
$\text{w}=(b,w_1,w_2,...,w_n{)}^T$，是列向量，表示回归系数。

---

总体思想

上述问题中如果利用线性回归对数据进行预测，但是对于类别仅有固定的两个的情况线性分类却无法正常使用，进一步思考，可以构造一种使得视力、身高以及体重成为线性关系的一组数据，利用线性回归求出来的$\text{y}$值对数据进行分类，那么现在最重要的就是找到这个$\text{y}$值。这种构造在理论上是可行的，将上面表格中的类别可以用具有线性关系的数据来代替，具体点可以说是年龄（18,19,20,10,12,13），那么很容易就可以通过线性回归的预测结果是否大于18岁来判定，但是这种构造很难找到。
退一步想，其实不用构造具体的数据(年龄)也可以解决，通过线性回归直接判定其结果是否大于18，而不用拘泥于某一组特定的数，基于这一思想，开始下面的理论探究。

---

推导过程

基于基本思想，利用线性回归公式${\text{w}}^T{\text{x}}_i$，由于0这个位置拥有较好的数学性质，不妨通过调节$\text{w}$中的b使其成为以0为分界点的回归函数(通过d=d-18，使以18岁为界改为以0为界)，当d≥0时为1类，反之为0类。
为了对数据进行训练，引入sigmoid函数。

介绍一下sigmoid函数，公式如下：
$$h(s)=\frac{1}{1+e^{-s}}$$
该函数有着良好的数学性质
1.$h(s)\in (0,1)$
2.${\lim }_{s\to +\infty }h(s)=1$,${\lim }_{s\to -\infty }h(s)=0$
3.$1-h(s)=h(-s)$
由性质1与性质2可以看出，该函数在任意点上可以表示某种概率，如果利用sigmoid构造概率函数则有
$$p(y|s)=\{\begin{array}{cc}h(s)& y=1\\ 1-h(s)& y=0\end{array}$$
等价于
$$p(y|s)=h(s{)}^y(1-h(s){)}^{1-y},y=1或0$$
接下来思考一下一种概率关系，当p(1|s)>p(0|s)时，那么可以说s为1类的概率大，比较有可能属于1类。由于h(s)以及1-h(s)都恒大于0，所以大于关系也可以转换为比例关系，即当$\frac{p(1|s)}{p(0|s)}>1$时，为1类，反之为0类。
$$\frac{p(1|s)}{p(0|s)}>1$$
$⟺$
$$\frac{h(s)}{1-h(s)}=\frac{1}{e^{-s}}=e^s>1$$
$⟺$
$$s>0$$
看到这步推导就可以明白这样构造概率函数的原因了，这种构造使概率关系与线性关系完美地连接起来，进一步地解释一下，令$s_i={\text{w}}^T{\text{x}}_i$，当$s_i>0$，即${\text{w}}^T{\text{x}}_i>0$时，那么一定有$p(1|s)>p(s)\Rightarrow h({\text{w}}^T{\text{x}}_i)>1-h({\text{w}}^T{\text{x}}_i)$,必定为1类，当${\text{w}}^T{\text{x}}_i>0$时，那么一定有$p(1|s)<p(s)\Rightarrow h({\text{w}}^T{\text{x}}_i)<1-h({\text{w}}^T{\text{x}}_i)$，必定为0类。此时概率关系（比较概率大小判断属于哪一类）与线性关系（总体思想中提到的${\text{w}}^T{\text{x}}_i>0$与${\text{w}}^T{\text{x}}_i<0$）表现出一致的结论。
另一个关键的部分是引入概率之后很方便我们使用极大似然估计法求解$\text{w}$的值,根据极大似然法有
$$P(\text{w})=\prod _{i=1}^{m}p(y_i,s_i)=\prod _{i=1}^{m}p(s_i)p(y_i|s_i)$$
假设$s_i$是等概率出现的(即$p(s_i),i=1,2,\cdots$的值全部相等)，那么原式可以写为如下形式，使$P(\text{w})$取最大值时不影响$\text{w}$的取值
$$P(\text{w};b)=\prod _{i=1}^{m}p(y_i|s_i)=\prod _{i=1}^{m}h(s_i{)}^{y_i}(1-h(s_i){)}^{1-y_i}$$
为方便求解，对该式取对数得到
$$L(\text{w})=\ln P(\text{w})=\sum _{i=1}^m{y}_{i}ln(\frac{h(s_i)}{1-h(s_i)})+ln(1-h(s_i))$$
我们的目标变为
$$\arg \underset{\text{w}}{\max }L(\text{w})$$
在求解过程中，我们倾向于求最小值
目标变为$$\arg \underset{\text{w}}{\min }-L(\text{w})$$

---

求解方法

$$\begin{array}{c}\begin{array}{rlr}L(\text{w})& =& \sum _{i=1}^m{y}_{i}ln(\frac{h(s_i)}{1-h(s_i)})+ln(1-h(s_i))\\ & =& \sum _{i=1}^m{y}_i{s}_i-\ln (e^{s_i}+1)\\ & =& \sum _{i=1}^m{y}_i{\text{w}}^T{\text{x}}_i-\ln (e^{{\text{w}}^T{\text{x}}_i}+1)\end{array}\end{array}$$
由于我们的求解目标是$-L(\text{w})$的最小值，为了方便表达，不妨令$L(\text{w})=-L(\text{w})$
那么
$$L(\text{w})=-L(\text{w})=\sum _{i=1}^m\ln (e^{{\text{w}}^T{\text{x}}_i}+1)-y_i{\text{w}}^T{\text{x}}_i$$
求梯度
$$▽L(\text{w})=\sum _{i=1}^m(h({\text{w}}^T{\text{x}}_i)-y_i){\text{x}}_i={\text{x}}^T({\text{y}}^{{}^{\prime }}-\text{y})\text{\hspace{1em}(1)}$$
$${\text{y}}^{{}^{\prime }}=(h({\text{w}}^T{x}_1),h({\text{w}}^T{x}_2),\cdots ,h({\text{w}}^T{x}_m){)}^T$$
求二阶梯度
$$▽{}^{2}L(\text{w})=\sum _{i=1}^m({\text{x}}_{i}h({\text{w}}^T{\text{x}}_i))h(-{\text{w}}^T{\text{x}}_i)){\text{x}}_i^T={\text{x}}^T\text{R}\text{x}\text{\hspace{1em}(2)}$$
其中$\text{R}$为mxm的对角方阵，对角元素${\text{R}}_{ii}=h({\text{w}}^T{\text{x}}_i)h(-{\text{w}}^T{\text{x}}_i)$。
牛顿法
牛顿法迭代公式
$${\text{w}}_{t+1}={\text{w}}_t-▽{{}^{2}L(\text{w})}^{-1}▽L(\text{w})$$
公式(1)与公式(2)带入求解
牛顿法在理论上可以求解，但是在应用中由于小数点位数的原因，往往会出现二阶梯度变为0矩阵的情况，导致无法继续迭代计算。

梯度下降法
梯度下降法迭代公式
$${\text{w}}_{t+1}={\text{w}}_t-\alpha ▽L(\text{w})$$
将公式(1)带入求解，得到$\text{w}$的值

import numpy as np
#simgoid函数
#直接return 1/(1+exp(-x))会报错
#由于x<0时,exp(-x)->无穷大，导致溢出
#当x<0时，按照exp(x)/(1+exp(x))计算，计算机可以超出位数的小数部分进行截断，使其继续保持计算 
def h(x):
	y=[]
	for ele in x:
		if ele>0:
			y.append(1/(1+np.exp(-ele)))
		else:
			y.append(np.exp(ele)/(1+np.exp(ele)))
	return np.array(y).reshape(x.shape)
#一阶梯度
def df1(x,y,w):
	return np.dot(x.T,h(np.dot(x,w))-y)
#二阶梯度
def df2(x,y,w):
	dr=h(np.dot(x,w))*h(-np.dot(x,w))
	m=dr.shape[0]
	r=np.zeros([m,m])
	for i in range(m):
		r[i,i]=dr[i]
	#print(np.dot(np.dot(x.T,r),x))
	return np.dot(np.dot(x.T,r),x)

#梯度下降法
def log_des(data,alp):
	dm,dn=data.shape
	w1=np.ones([dn,1])
	w0=np.zeros([dn,1])
	y=data[:,-1:]
	x=np.append(np.ones([dm,1]),data[:,:-1],axis=1)
	while sum(abs(w1-w0))[0]>0.00001:
		w0=w1
		w1=w0-alp*df1(x,y,w0)
	return w1

#牛顿法
def log_newton(data):
	dm,dn=data.shape
	w1=np.ones([dn,1])
	w0=np.zeros([dn,1])
	y=data[:,-1:]
	x=np.append(np.ones([dm,1]),data[:,:-1],axis=1)
	while sum(abs(w1-w0))[0]>0.0001:
		w0=w1
		#由于运算位数原因，在迭代过程中二阶梯度变成了0矩阵，导致逆不存在，无法计算
		w1=w0-np.linalg.solve(df2(x,y,w0),df1(x,y,w0))
	return w1
data=np.array([[3.8,180,85,1],[3.7,189,81,1],[3.8,190,100,1],[5.2,160,50,0],[5.1,166,70,0],[5.0,170,60,0]])
#提取前5行作为训练集
w=log_des(data[:-1,:],0.01)
print("梯度下降法参数w结果：\n",w)
#将所有数据预测一下
print("预测结果：\n")
for ele in np.dot(np.append(np.ones([data.shape[0],1]),data[:,:-1],axis=1),w):
	if ele>0:
		print(1)
	else:
		print(0)
#print("牛顿法结果：",log_newton(data))

利用sklearn做一下，没找到输出训练参数的函数

import numpy as np
from sklearn.linear_model import LogisticRegression
tdata=np.array([[3.8,180,85,1],[3.7,189,81,1],[3.8,190,100,1],[5.2,160,50,0],[5.1,166,70,0],[5.0,170,60,0]])
model=LogisticRegression()
#对数据进行训练
model.fit(tdata[:-1,:-1],tdata[:-1,-1:].ravel())
#预测结果
print(model.predict(tdata[:,:-1]))
#没找到如何输出训练参数