linear regression and logistic r

2018-08-23 本文已影响31人冒绿光的盒子

①linear regression

target function的推导

线性回归是一种做拟合的算法：

通过工资和年龄预测额度，这样就可以做拟合来预测了。有两个特征，那么就要求有两个参数了，设置

一般误差不会太多，都是在0上下浮动的，所以0附近是最高的概率。
使用linear regression要满足三个条件：
1.独立，每一个样本点之间都要相互独立。
2.同分布，他们的银行是一样的，使用的算法是一样的。
3.误差都服从高斯分布。

由于误差是服从高斯分布的，自然有： $P(\theta_0)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{\theta_0^2}{2\sigma^2})$
上述式子有： $y = \theta^Tx+\theta_0$
$P(\theta_0)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y_i-\theta^Tx_i)^2}{2\sigma^2})$
需要求参数 $\theta$ ，自然就是最大似然函数：
$L(\theta)=\prod_{i=1}^m\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y_i-\theta^Tx_i)^2}{2\sigma^2})$
log化简：
$logL(\theta) = log\prod_{i=1}^m\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y_i-\theta^Tx_i)^2}{2\sigma^2})$
$=mlog\frac{1}{\sqrt{2\pi}\sigma}+\sum_{i=1}^m(-\frac{(y_i-\theta^Tx_i)^2}{2\sigma^2})$
$=mlog\frac{1}{\sqrt{2\pi}\sigma}-\frac{1}{2\sigma^2}\sum_{i=1}^m(y_i-\theta^Tx_i)^2$
最大化这个似然函数那就是最小化 $\frac{1}{2\sigma^2}\sum_{i=1}^m(y_i-\theta^Tx_i)^2$ 。因为前面都是常数，可以忽略。而这个式子稍微变一下就是最小二乘法了。
最小二乘法： $f = \frac{1}{2}\sum_{i=1}^m(h_{\theta}(x_i)-y_i)^2$ ，不用上面化简出来的式子非要加上 $\frac{1}{2}$ 是因为求导之后2是可以和0.5抵消的。

目标函数的化简

化简一下：
$J(\theta)=\frac{1}{2}\sum_{i=1}^m(h_{\theta}(x_i)-y_i)^2$
$=\frac{1}{2}(x\theta-y)^T(x\theta-y)$
求极小值，直接就是导数为0即可，所以对 $\theta$ 求导数即可。
$\bigtriangledown_{\theta}J(\theta) = \bigtriangledown_{\theta}(\frac{1}{2}(\theta^Tx^T-y^T)(x\theta-y))$
$= \bigtriangledown_{\theta}\frac{1}{2}(\theta^Tx^Tx\theta-\theta^Tx^Ty-y^Tx\theta+y^Ty)$
$=x^Tx\theta-x^Ty$
直接令等于0： $\theta = (x^Tx)^{-1}x^Ty$ 。
这样就是得到了最终的目标函数，一步就可以到达最终的optimal value。这就是linear regression。所以线性回归是有解析解的，直接一条公式就可以求出来的。注意上面的公式 $J(\theta)=\frac{1}{2}\sum_{i=1}^m(h_{\theta}(x_i)-y_i)^2$ 的 $\theta^Tx$ 是不包含bias偏置值，偏置项就是误差，代入高斯函数推导出来的。

②logistic regression

target function 的推导

首先要提一个函数，sigmoid函数： $f(x) = \frac{1}{1+e^{-z}}$
这个函数之前被用来做神经网络的激活函数，但是它有一个缺点。

离原点处越远梯度就越趋近于0，在神经网络里面就会出现一个梯度消失的问题，所以后来就用relu或者tanh函数代替了。这个函数的值域是在

所以梯度更新公式：

有三个不同的sigmoid函数，两个式子解不了，因为sigmoid函数不是线性的，如果是线性的，那么可以吧

\theta

提出来，这样就有解了。其实如果sigmoid去掉，仅仅把

sigmoid(\theta^Tx) = \theta^Tx

那么就和linear regression一样了。所以是没有解析解的，主要的原因就是因为sigmoid函数是一个非线性的函数。

当标签不同的时候另一种函数形式

计算分数函数一样的： $score = W^Tx$
之前了解的preceptron：

直接把score分数通过一个sign函数即可。0/1错误。
linear regression:

线性回归就是去掉了sign函数，使其成为一个线性函数，error function = square
logistic regression：

常规操作，就是要找到error function先。按照刚刚的函数分布可以得到：

极大似然函数：

sigmoid函数有一个比较牛逼的性质。 $1-\frac{1}{1+e^{-x}}=\frac{e^{-x}}{1+e^{-x}}=\frac{1}{\frac{1}{e^{-x}}+1} = \frac{1}{1+e^x}$

最大值一般比较难求，加上负号变成一个最小化的问题，化简一下：

其实这个就是交叉熵函数，和上面推导出的：

是一个东西，整合起来了而已，负号提前，所以形式有点不一样。最小化error function：

gradient decent，和刚刚的方法是一样的。

上面的 $\theta(yw^Tx)$ 可以看做是-yx的线性加权，要求就是加权和之后为0，而如果是线性可分的话，则只要 $\theta$ 为0就可以了；根据sigmoid函数的特性，为0就相当于是 $-yw^Tx$ 要 << 0，即 $-yw^Tx$ >> 0，这就尴尬了，需要保证全部的yx都同号，都是线性可分的，这样其实是很难做到的，所以我们转换一个思路，用梯度下降解决。归根到底，还是sigmoid函数的非线性。

Stochastic Gradient Descent

随机梯度下降，以往的经验是全部一起做一次更新，如果计算量非常大，那么计算复杂度很高，另外如果是做online learning的时候也不方便，因为这个时候数据不是一个betch的过来了，而是几个几个的了。

优点：简单，如果迭代次数够多的话是可以比得上average true gradient的效果的。
缺点：迭代的方向会不稳定，可能会左拐右拐的。average true gradient是查找的当前一阶导最适合的方向走的，而SGD是直接随机一个方向走。

他和PLA其实很像，而在理解上面其实也是类似，他们都选取一个点进行调整。而SGD ≈ soft PLA，PLA是错了之后才纠正，而SGD会看下你错了多少，错了越多，意味着你的wx越大，所以纠正的越多，是动态可变的，所以也叫做soft PLA。

③线性模型error function的对比

三个比较简单算法：PLA，linear regression，logistic regression。他们勇于分类的时候：

square function对于分类来说其实不太合理的，分类正确了，应该越远越好才对，但是square function是越远错误就越大，是不合理的，logistics就更合理了，错误的越错就越大正确的就小，所以linear regression适合回归而不是分类。
可以看到ce和err0/1是有焦点的，我们乘上一个数字使他相切：

根据VC bound的理论：

所有logistic regression是可以作为分类的，而且他的分类效果要比linear regression好，首先直观理解错误，他比linear regression更合理，其次，他的VC bound比linear regression要小，这就证明了Ein ≈ Eout的概率会更高，泛化能力更强。

④Nonlinear Transformation

对于线性可分的情况来说，几乎是不用对x做什么预处理就可以直接使用模型进行分类，但是如果是对于非线性的模型，上面的方法就有点吃力了，他们都是线性分类，直接在model上改进有点困难，所以在数据上进行处理。

这种是绝逼找不到一条直线分开的，要完成这种任务就需要找一个非线性的模型分开，比如看起来这就像是一个圆：

把上面的式子变成我们认识的形式：

支持向量机也用到了这种思想，核函数就是映射到高维度空间里面然后进行切片操作。
Nolinear Transformation的方法不止一个：

目前我们所讨论的都是过原点的，如果是不过原点的话，那么他们的VC dimension就会增加。比如这种

整个过程就是通过映射关系，换个空间去做线性分类，重点包括两个：
特征转换
训练线性模型

Price of Nonlinear Transform

首先，用十二指肠想都知道： $d_{vc}^{nonlinear} >=d_{vc}^{linear}$
如果是d维的x做二次的扩展，那么有

求上界，如果阶数更高，假如阶数为Q，对于x是d维，那么对于z空间就是

这种特征变换会导致计算空间变大，计算复杂度变大。另外，可以看到随着Q增大，W也会增大，这样就导致了VC dimension会增大，而W的秩其实就是VC维，所以，如果Q越大，就会导致泛化能力变差。

接下来讨论一下x到z多项式的变换：
一维：

对应上图：

所以一目了然，代价就是 $E_{out},E_{in}$ 之间的差距会越来越大， $E_{in}$ 最后就不能再代表 $E_{out}$ 了。随着变换多项式的阶数增大，虽然逐渐减小，但是model complexity会逐渐增大，造成很大，所以阶数不能太高。如果选择的阶数很大，确实能使接近于0，但是泛化能力通常很差，我们把这种情况叫做tempting sin。所以，一般最合适的做法是先从低阶开始，如先选择一阶hypothesis，看看是否很小，如果足够小的话就选择一阶，如果大的话，再逐渐增加阶数，直到满足要求为止。也就是说，尽量选择低阶的hypothes，这样才能得到较强的泛化能力。模型复杂度越高出现的泛化能力问题，就是过拟合，在训练数据上表现很好，Ein = 0，但是在测试数据上就很差，Ein << Eout。

⑤Overfitting and Regularization

overfiting，过拟合，就是上面所描诉的情况了， $E_{in}$ 和 $E_{out}$ 相差太远了，局部已经不能再代表全局了，主要原因就是模型复杂度太大了，VC bound太大，限制不了 $E_{out}$ 。overfitting就是训练过程中翻车了，原因：①太快了，模型复杂度太大。②近视，看的路太短，样本少。③这路不行，弯弯曲曲的，noise太多，噪音太大。
对应的解决方法其实很多，但主要就是regularization了：

虽然fitting出来的结果完全符合了数据点，但是模型本身没有这么复杂，所以产生了过拟合。
既然模型太复杂了，那么简化一下模型：

所以过拟合了，我们可以简化复杂度，复杂度的代表系数其实就是VC dimension，也就是W的个数。比如一个十阶的hypothesis：

这样就把W限定在了一个以根号C的一个球形里面，两个W才是一个圆，多个w就是立体的球了。

可以看到梯度的方向是差不多对准了 $W_{lin}$ （最好的W），只要在切线方向有分量，那么就会沿着切线的方向滑动。但是可以看到，这梯度的方向只是差不多对准了 $W_{lin}$ 而不是准确的对齐，因为是直接的求导，你不可以保证这个式子就是一个凸优化问题，它可能是多个凸优化多个山谷的，可能只是看到了比较近的一个山谷，但是大致方向是正确的。所以如果梯度没有和分量w平行，那么会一直移动， $w,normal$ 可以看做是 $w,w^T$ ，梯度就是 $-△E_{in}$ 负方向嘛，当平行时，按照图中的推理： $△E_{in}+\frac{2\lambda}{N}w = 0$ ， $△E_{in}$ 就是square error function求导，所以有 $\frac{2}{N}(z^Tzw-z^Ty)+\frac{2\lambda}{N}w=0$ ，最后就推导出 $w = (z^Tz+\lambda I)^{-1}z^Ty$
这种线性回归被看成是origin linear regression的进阶版，也叫ridge regression。
上面的公式也可以直接用拉格朗日乘子法推出来：

可以看到λ越大w越小，对于模型的惩罚也就越大，所以也叫做正则化惩罚项。这种惩罚会是的曲线越来越平滑，后面还会讲到另一种regularization。

对于regularization的VC 理论解释

正则化出来过后的linear regression就是ridge regression，根据之前的VC bound理论： $E_{out}(w) <= E_{in}(w)+Ω(H)$ ，其中 $Ω(H)$ 是复杂度， $Ω(w)=w^Tw$ 代表的是单个hypothesis的复杂度，而 $Ω(H)$ 代表的整个hypothesis set的复杂度，所以 $Ω(w)$ 是被包含在 $Ω(H)$ 里面的，所以相对来说 $E_{in}(w)+w^Tw < E_{out}(w)+Ω(H)$ ，相对来说会和 $E_{out}$ 更加接近。

而对于VC dimension，既然w被限制了，那么ridge的肯定比origin的小了。

General Regularizers

两个比较常用的：
ridge：L2范数， $Ω(w) = |w|^2$
lasso：L1范数， $Ω(w) = |w|$
ridge推导过了，来看看lasso的：

他的正常方向应该是垂直于边界的，红色的sign就是方向，而对于在边界上，只要sign和△Ein不在一条直线上，那么在正方形的边界上就一定有一个类似于刚刚ridge绿色分量的分量存在，如上图，就会向上方移动，而在边角点其实是不可微分的，所以一般结果会聚集在边角的周围，所以得到的是一个稀疏矩阵，一些是0，一些不是零。结果会聚集在顶点周围，优点就是计算很快了。
上图是用直观的理解来解释regularization，前面的是用VC demension来解释，再用数值分析的角度看一下：
比如：方程式①：
5x + 7y = 0.7
7x + 10y = 1
x = 0 , y = 0.1
方程式②：
5x + 7y = 0.69
7x + 10y = 1.01
x = -0.17 y = 0.22
只要有一点微小的扰动，结果就会发生很大的变化，这就是病态矩阵，如果 $X^TX$ 是病态矩阵，那么两次的结果可能会相差很大。是不是病态矩阵可以用条件数来判断。用条件数来描述病态矩阵其实还是太抽象了。奇异矩阵就是不存在逆矩阵的方阵。什么样的方阵不存在逆矩阵？首先行列式为0，自然就牵扯出来非满秩，特征值只和为0。而对于近似奇异矩阵，他的行列式很接近于0，所以 $\frac{1}{|A|}$ 是一个很大的数字，根据求逆矩阵的公式，逆矩阵是可以通过伴随矩阵和行列式求的，所以自然差别就很大了，所以病态矩阵也就是近似于奇异矩阵的矩阵了。那么知道原因了，我们要做的就是远离奇异矩阵。而事实上，ridge regression的结果：

加上一个对角线的值其实就是远离奇异矩阵。
最后再看一个图：

lasso可以使得权值衰减到0，很容易可以和顶点相切，所以会使得权值衰减导0；但是ridge只会减少，而不会衰减导0，所以lasso是可以用来控制权值的数量的。lasso同上也是满足了拉普拉斯分布，是一个稀疏矩阵，而ridge是高斯分布。

summary

加入正则项是为了避免过拟合,或解进行某种约束,需要解保持某种特性
①L1正则可以保证模型的稀疏性，也就是某些参数等于0,L1正则化是L0正则化的最优凸近似，比L0容易求解，并且也可以实现稀疏的效果,
②L2正则可以保证模型的稳定性，也就是参数的值不会太大或太小.L2范数是各参数的平方和再求平方根，我们让L2范数的正则项最小，可以使W的每个元素都很小，都接近于0。但与L1范数不一样的是，它不会是每个元素为0，而只是接近于0。越小的参数说明模型越简单，越简单的模型越不容易产生过拟合现象。
③在实际使用中，如果特征是高维稀疏的，则使用L1正则；如果特征是低维稠密的，则使用L2正则。
④L2不能控制feature的“个数”，但是能防止模型overfit到某个feature上；相反L1是控制feature“个数”的，并且鼓励模型在少量几个feature上有较大的权重。
总结一下之前学过的，第十八篇博客VC dimension是机器学习的理论保证：

霍夫丁不等式保证了一个hypothesis发生坏事的概率是很小的，Multi霍夫丁不等式保证了这个hypothesis set对于所有的数据集发生的坏事的概率满足的不等式，VC dimension则保证了这个hypothesis set里面发生的坏事的概率是很小的，这就保证了泛化能力：Ein ≈ Eout。
之后介绍了三种线性模型（PLA比较简单，不会再讲了）：

而对于这三种模型，我们又给出了三种工具：

Feature Transform：通过低维不可线性可分转换到高维线性可分来解决Nonlinear的情况。
Regularization：通过增加正则化惩罚项来解决过拟合的问题。
Validation：用于对模型参数的选择，第一次选择：是选择上面模型。第二次选择就是对模型参数的选择了，这就用Validation来帅选。
最后是三个比较常用的锦囊妙计：

⑥logistics regression代码实现

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

def predict(w, x):
    h = x * w
    h = sigmoid(h)
    if h > 0.5:
        return int(1)
    else:
        return int(0)

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

def error_rate(h, label):
    m = np.shape(h)[0]
    sum_err = 0.0
    for i in range(m):
        if h[i, 0] > 0 and (1-h[i, 0]) > 0:
            sum_err += (label[i, 0] * np.log(h[i, 0]) + (1-label[i, 0])*np.log(1-h[i, 0]))
        else:
            sum_err += 0
    return sum_err/m

def lr_train_bgd(feature, label, maxCycle, alpha, df):
    n = np.shape(feature)[1]
    m = np.shape(feature)[0]
    w = np.mat(np.random.rand(n,1))
    i = 0
    while True:
        i += 1
        h = sigmoid(feature * w)
        err = label - h
        if i % 100== 0:
            print('error : ', error_rate(h, label))
            d = 0
            scores = []
            for i in range(m):
                score = predict(w, feature[i])
                scores.append(score)
                if score == label[i]:
                    d += 1
            print('train accuracy : ', (d/m)*100, '%')
            if (d/m)*100 >= 90:
                for i in range(m):
                    if df.iloc[i, 2] == 1:
                        c = 'red'
                    else:
                        c = 'blue'
                    plt.scatter(df.iloc[i, 0], df.iloc[i, 1], c=c)
                x = [i for i in range(0, 10)]
                plt.plot(np.mat(x).T, np.mat((-w[0]-w[1]*x)/w[2]).T, c = 'blue')
                plt.show()
                return

        w += alpha * feature.T * err

def loadData(filename):
    df = pd.read_csv(filename, sep='    ', names=['a', 'b', 'label'])
    n, m = df.shape
    features = []
    labels = []
    for i in range(n):
        feature = []
        feature.append(int(1))
        for j in range(1,m):
            feature.append(df.iloc[i, j])
        if df.iloc[i, m-1] == -1:
            labels.append(0)
        else:
            labels.append(1)
        features.append(feature)

    for i in range(n):
        if df.iloc[i, 2] == 1:
            c = 'red'
        else:
            c = 'blue'
        plt.scatter(df.iloc[i, 0], df.iloc[i, 1], c = c)
    plt.show()
    return np.mat(features), np.mat(labels).T, df

if __name__ == '__main__':
    f, t, df = loadData('../Data/testSet.txt')
    lr_train_bgd(f, t, 10000, 0.001, df)