机器学习->监督学习->logistic回归,softMax回归

本篇博文来总结一下回归模型里面两个非常重要的模型。

• logistic回归
• softMAX回归

Logistic回归

logistics回归虽然有“回归”两字但是却是分类模型，并且是个二分类模型。

logistics回归是个线性分类模型，有着线性的决策边界，但是有着非线性的激活函数去估计后验概率。

下面就从数学层面讲讲logistics回归。
首先介绍下sigmoid函数

其函数图像：

显然sigmoid函数只能取得0到1之间的数。
线性分类有着线性的决策边界，类似于下图：

对于线性边界，边界决策函数：

构造预测函数：

现在以此作为二分类分类标准，下面分别是为类别1，类别0的概率：

综上：

在衡量模型好坏方面除了构造合理的损失函数外还可以通过最大似然估计，即是用似然函数来进行参数估计。

似然函数（独立作积）：

最大似然函数（目标函数）：

上式其实就是将$cross-entropy$ 的公式取负：

$$L =- \sum p(x)log(q(x))$$
用来衡量两种分布的拟合程度，如果 $p(x)、q(x)$ 为两种相同的分布，则其 $cross\ entropy$ 为0。

最大化$MLE\ loss$ 其实就是最小化$cross\ entropy$。

我们这里是要求出似然函数最大时对应的参数，一般来说对似然函数求导，其导函数为0时对应的参数就是我们要求的参数，但是这里很难对似然函数直接求导使其为0。

同线性回归一样，我们可以在目标函数后面加上正则项来控制模型的复杂程度，来防止过拟合，可以加上L1或L2正则项，具体可以看这篇博文机器学习->监督学习->线性回归（LASSO,Ridge,SGD） 里面相关内容。

最大似然函数对求导，然后通过梯度下降法来求梯度方向：

这里需要注意：当我们以损失函数进行参数估计时（例如线性回归里面的最小二乘估计），因为是求损失函数最小时的参数（是一个不断下降地逼近最低值），故采用梯度下降法，当我们用最大似然函数进行参数估计时，是求似然函数最大时对应的参数（是一个不断上升地逼近最高值），那么这时时梯度上升法。其实道理都是一样，只是前面正负号的问题。

我们发现通过最大似然估计得出的梯度方向还是一种的形式，跟之前线性回归得出的梯度公式一模一样，那么这是为什么呢？这里需要介绍下一个重要概念 指数族分布，有关这个概念更详细的说明请看这篇博文指数族分布。那么哪些分布属于指数族分布呢？

• 高斯分布属于指数族分布
• 二项分布属于指数族分布
• 泊松分布属于指数族分布

只要样本的分布属于指数族分布，那么在用最大似然估计求梯度方向时，最后的形式都是一个，并且得出的模型都是线性模型，只不过是一种广义线性模型（GLM）。

我们再来明确一下什么叫广义线性模型（GLM）

• y不再只是正态分布，而是扩大为指数族中 的任一分布
• 变量
  连接函数g
• 连接函数g单调可导
  如Logistic回归中

在梯度下降法中每次迭代都要代入全部数据进行更新参数，在批梯度下降法中，每次迭代随机选择一部分数据样本来更新参数：

当数据量小时还好，当数据量很大时，这样迭代方法耗时耗力。
由此产生mini_batch梯度下降法（mini_batchGD），就是每次迭代随机选择训练集的中一个数据样本或者一批数据样本代入更新参数。

随机梯度法虽然大大减少了迭代所需的时间，但是偶尔却会使预测结果并不十分准确。

def logistic_regression(data,alpha,lamda):
    n=len(data)-1　## ｎ个特征
    w=[0 for x in range(n)]
    w1=[0 for x in range(n)]
    for times in range(10000):
        for d in data:
            for i in range(n):
                w1[i]=w[i]+alpha*(d[n]-h(w,d))*d[i]+lamda*w[i] ## 加上了一个正则项
            for i in range(n):
                w[i]=w1[i]
    return w

Logistic Regression + Square Error

首先，线性回归中其$\hat y$ 为连续值，而逻辑回归中$\hat y$ 为离散的数值。

我们在线性回归中，讲到用$square\ Error$ 作为损失函数，那么为什么在逻辑回归中，不采用$square\ loss$ 呢？
没事，我们先尝试采用$square\ loss$ 应用到逻辑回归中，则有：

由上面的推导可以看出：
当$\hat y==0$ 时：

• 计算出的$f(x) == 1$ 时，说明距离最优点还很远，但是此时计算出的导数为0，参数几乎不更新了。
• 计算出$f(x) == 0$ 时，说明已经在最优点了，此时导数为0，看起来较为合理。

$\hat y == 1$ 时，同样有这种缺陷。

所以在逻辑回归中，不太适用$square\ loss$。

由上面这张图来看，$cross-entropy$ 在距离最优点较远的地方时，斜率较大，参数更新较快；而$square-Error$ 在距离最优点较远地方，变化比较平滑，更新较慢。

Logistic回归的损失函数

在逻辑回归中，，由此可得：

我们要取这个对数似然函数的最大值。那么（negative Log Likelihood，简称NLL）就可以认为是逻辑回归的损失函数。

因为这个最后的形式过于复杂，我们可以做个简单的变换使得
很明显当时，，。

Limitation of Logistic Regression

逻辑回归毕竟是线性的分类，还是存在天然无法的一些问题，例如：

显然对于上面，对于上面这类$DataSet$，逻辑回归不好解决。那么如果非要用逻辑回归解决，该如何是好呢？
先做$Feature\ transformation$，将某一类的样本分到直线的一侧，例如：

SoftMax回归

对Logistic回归有了个深刻的认识后，softMax回归就简单多了，只不过由二分类变成了多分类。

def soft_max(data,K,alpha,lamda):
    n=len(data)-1 ##　样本特征维度
    w=np.zeros((K,n)) ##当前权重矩阵：每个类别有各自的权重向量
    w1=np.zeros((K,n)) ## 临时权重矩阵，用于存放临时更新的权重
    for time in range(1000):
        for d in data:
            x=d[:-1] ## 输入特征向量
            for k in range(K):
                y=0  ## 默认情况下，这个样本属于ｋ类的期望（概率）为０
                if int(d[-1]+0.5)==k:
                    y=1 ## 在0.5误差内，在认为该样本属于ｋ类的期望为１
                p=predict(w,x,k) ## 输出预测该样本为k类的期望
                g=(y-p)*x ## 梯度
                w1[k]=w[k]+alpha*g+lamda*w[k]
            w=w1.copy()
    return w

根据以上所学做一道题目：

题目和数据链接

#coding:utf-8
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.patches as mpatches


def normalization(X):
    '''
    X标准化：以一列作为一个单位，标准化公式(x-该列均值）/（该列标准差）
    因为给的X有两列，Y一列，则对应的theta有三列，其中theta0对应的X默认是1，故从X的第二列开始标准化。
    '''
    for i in range(1,int(X.shape[1])):
        std,mean=np.std(X[:,i],ddof=1),np.mean(X[:i])##标准差公式里面默认(ddof=0)除以N，就是数组长度,ddof=1就是除以(N-1)
        X[:,i]=(X[:,i]-mean)/(std+1e-10)
    return X

def sigmoid(X,T):
    '''
    :param X: 特征矩阵
    :param T: 参数矩阵
    :return: 返回预测概率
    '''
    return 1.0/(1+np.exp(-(X.dot(T))))

def loadData():
    data = pd.read_csv('iris.data', header=None)
    data[4] = pd.Categorical(data[4]).codes
    data=data[data[4]!=2]
    x, y = np.split(data.values, (4,), axis=1)
    x = x[:, :2]
    ones=np.ones(shape=(x.shape[0],1))
    x=np.c_[ones,x]##在第一列插入全1
    x=normalization(x)## 数据归一化，这一步非常重要
    return x,y

def MSError(y,y_hat):
    '''
    :param y:真实y值
    :param y_hat: 预测ｙ值
    :return: 返回均方误差
    '''
    return ((Y-y_hat).transpose().dot((Y-y_hat)))/(2*y.shape[0])

def GD(X,Y):
    epsilon=1e-8
    alph=0.001
    error1=0
    cnt=0
    m=np.shape(Y)[0]
    T=np.ones((X.shape[1],1))
    while True:
        cnt=cnt+1
        y_hat=sigmoid(X,T)
        error=y_hat-Y
        T=T-alph*X.transpose().dot(error)## 根据error×feature形式更新参数矩阵
        error0=error1
        error1=MSError(Y,y_hat)##计算均方误差
        if abs(error1-error0)<epsilon:##如果上一次的均方误差和当前次的均方误差很小（可视为误差收敛），则停止
            break
        plot2D(X,Y,T,cnt,error1)
    print "iterations",cnt
    return T

def plot2D(X,Y,T,cnt,error1):
    axarr[0].cla()
    m=np.shape(Y)[0]
    cm_dark = mpl.colors.ListedColormap(['b', 'g'])
    axarr[0].scatter(X[:,1],X[:,2],c=Y,cmap=cm_dark)
    x1_min, x1_max = X[:, 1].min()-3, X[:, 1].max()+3   # 第0列的范围
    x1=range(int(x1_min),int(x1_max),1)
    T=np.array(T.transpose())[0]
    x2=[-(T[0]+i*T[1])/T[2] for i in x1]#画样本分割线，分割线上对应的Y值为0；t0+t1*x1+t2*x2=0
    axarr[0].plot(x1,x2)
    axarr[0].set_xlim([x1_min, x1_max])
    axarr[0].set_ylim([-5, 5])
    axarr[0].set_xlabel('x1')
    axarr[0].set_ylabel('x2')
    axarr[0].set_title('data')
    axarr[1].plot(cnt, error1, 'or')
    axarr[1].set_xlim([0, 200])
    axarr[1].set_ylim([0, 0.8])
    axarr[1].set_xlabel('step')
    axarr[1].set_ylabel('cost')
    axarr[1].set_title('cost')
    fig.canvas.draw()


if __name__=='__main__':
    fig, axarr = plt.subplots(1, 2, figsize=(15, 5))
    fig.show()
    X, Y = loadData()
    T = GD(X, Y)