多元Logistic回归模型 多元logistic回归模型拟合信息

有监督模型

1. 线性回归

1.1 线性模型

xi和h0是已知的观测值，需要求得theta的参数值，这一步其实就是在拟合一个截断面。当求得theta后，只需要给定xi的值，就可以推断出h0的值，起到了预测的作用。

1.2 如何求得theta

1.2.1 似然函数

如果theta是被已经求得，那么预测如下公式所示

。y_pred会与真实值有一定的偏差，我们将偏差记作ε，此时真实值与预测值之间的关系可以表示为（即似然函数）：

。

epsilon表示误差，当误差变小时，p(yi|xi;theta)的值在变大，那么我希望在已知很多样本的x和y值时总体误差越小越好，说明拟合接近真值，那么对应到p(yi|xi;theta)就是越大越好，又因为各样本之间独立，总体误差小=总体p(yi|xi;theta)大=

连乘值大。

概率分布函数与概率密度函数可以参考如何理解概率分布函数和概率密度函数？_mjiansun的专栏-CSDN博客_概率分布与概率密度

对上述几个公式的转变进行说明：

公式是根据高斯概率密度函数所得，高斯概率密度函数公式为：

其中μ表示均值为0，σ表示标准差为固定值，带入偏差ε后即可得到公式。

上述公式可以看成已知theta和x值，关于真实值y的分布函数。那么此时的概率密度函数就是上式。

上面L(theta)能直接看成y1y2等同时出现的情况即L(y1y2...)，又由于独立同分布L(y1y2...)=L(y1)L(y2)L(y3)L(y4)...这就有了下面的公式。由于表达式中的x和y是观测值已经给定，theta未知是要求的值，那么此时就可以写作关于theta的联合概率密度函数。

1.2.2 对数似然与最小二乘

由于L(theta)是乘积的形式，为了方便计算，常规处理加对数，即得到对数似然。

其中公式中的第一项除了m都是定值，可以忽略。第二项的σ也是定值，可以忽略，其实1/2也是定值本质可以忽略，但是为了方便后续计算保留了1/2。将保留的部分单独提取出来，就形成了最小二乘的形式。

对于对数似然函数其值越大越好，其实就是求

的最大值，去掉前面的负号，那么就是求

的最小值，从而就有了最小二乘的形式。

最小二乘

此时的目标为最小化上述公式。

1.2.3 theta解析式的求解过程（第一种方式）

 

梯度具体如何求，请参考矩阵求导术（上）_mjiansun的专栏-CSDN博客

让梯度公式等于0，可求得驻点，即可认为此结果为参数最优解，到这里关于theta的解析式已经求出来。

1.2.3.1 针对解析式几个常见的概念和疑问

（1）线性模型怎么得到曲线？

线性是指a0，a1，a2等构成的参数空间与最终的求和运算是线性的。

1阶线性模型：y=a1*x + a0*1

X = (x1, 1),(x2, 1),(x3, 1),(x4, 1)....

Y = y1,y2,y3,y4,....

2阶线性模型：y=a2*x^2 + a1*x + a0*1

X = (x1^2, x1, 1),(x2^2, x2, 1),(x3^2, x3, 1),(x4^2, x4, 1)....

Y = y1,y2,y3,y4,....

3阶线性模型：y=a3*x^3 + a2*x^2 + a1*x + a0*1

X = (x1^3, x1^2, x1, 1),(x2^3, x2^2, x2, 1),(x3^3, x3^2, x3, 1),(x4^3, x4^2, x4, 1),....

Y = y1,y2,y3,y4,....

上述所谓的1阶，2阶，3阶的本质只是x的变化，可以将x^2，x^3等看成x2,x3即可，将其看成一种新的特征，那么依旧对应了

该公式的计算，只不过这里的x1，x2，x3等可以通过x指数次方求出来而已。 

利用上述

模型和一系列x和y就能求得theta值，也变向求出了模型。注：这里的theta的数量与特征的数量相等。

代码

#!/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.exceptions import ConvergenceWarning
import matplotlib as mpl
import matplotlib.pyplot as plt
import warnings
# import seaborn


def xss(y, y_hat):
    y = y.ravel()
    y_hat = y_hat.ravel()
    # Version 1
    tss = ((y - np.average(y)) ** 2).sum()
    rss = ((y_hat - y) ** 2).sum()
    ess = ((y_hat - np.average(y)) ** 2).sum()
    r2 = 1 - rss / tss
    # print 'RSS:', rss, '\t ESS:', ess
    # print 'TSS:', tss, 'RSS + ESS = ', rss + ess
    tss_list.append(tss)
    rss_list.append(rss)
    ess_list.append(ess)
    ess_rss_list.append(rss + ess)
    # Version 2
    # tss = np.var(y)
    # rss = np.average((y_hat - y) ** 2)
    # r2 = 1 - rss / tss
    corr_coef = np.corrcoef(y, y_hat)[0, 1]
    return r2, corr_coef


if __name__ == "__main__":
    warnings.filterwarnings(action='ignore', category=ConvergenceWarning)
    np.random.seed(0)
    np.set_printoptions(linewidth=300, suppress=True)
    N = 9
    x = np.linspace(0, 6, N) + np.random.randn(N)
    x = np.sort(x)
    y = x**2 - 4*x - 3 + np.random.randn(N)
    x.shape = -1, 1
    y.shape = -1, 1

    models = [Pipeline([
        ('poly', PolynomialFeatures()),
        ('linear', LinearRegression(fit_intercept=False))]),
        Pipeline([
            ('poly', PolynomialFeatures()),
            ('linear', RidgeCV(alphas=np.logspace(-3, 2, 10), fit_intercept=False))]),
        Pipeline([
            ('poly', PolynomialFeatures()),
            ('linear', LassoCV(alphas=np.logspace(-3, 2, 10), fit_intercept=False))]),
        Pipeline([
            ('poly', PolynomialFeatures()),
            ('linear', ElasticNetCV(alphas=np.logspace(-3, 2, 10), l1_ratio=[.1, .5, .7, .9, .95, .99, 1],
                                    fit_intercept=False))])
    ]
    mpl.rcParams['font.sans-serif'] = ['simHei']
    mpl.rcParams['axes.unicode_minus'] = False

    plt.figure(figsize=(15, 10), facecolor='w')
    d_pool = np.arange(1, N, 1)  # 阶
    m = d_pool.size
    clrs = []  # 颜色
    for c in np.linspace(16711680, 255, m, dtype=int):
        clrs.append('#%06x' % c)
    line_width = np.linspace(5, 2, m) * 0.7
    titles = '线性回归', 'Ridge回归', 'LASSO', 'ElasticNet'
    tss_list = []
    rss_list = []
    ess_list = []
    ess_rss_list = []
    for t in range(4):
        model = models[t]
        plt.subplot(2, 2, t+1)
        plt.plot(x, y, 'ro', markersize=7, zorder=N, mec='k')
        for i, d in enumerate(d_pool):
            model.set_params(poly__degree=d)
            model.fit(x, y.ravel())
            lin = model.get_params('linear')['linear']
            output = '%s：%d阶，系数为：' % (titles[t], d)
            if hasattr(lin, 'alpha_'):
                idx = output.find('系数')
                output = output[:idx] + ('alpha=%.6f，' % lin.alpha_) + output[idx:]
            if hasattr(lin, 'l1_ratio_'):   # 根据交叉验证结果，从输入l1_ratio(list)中选择的最优l1_ratio_(float)
                idx = output.find('系数')
                output = output[:idx] + ('l1_ratio=%.6f，' % lin.l1_ratio_) + output[idx:]
            print(output, lin.coef_.ravel())
            x_hat = np.linspace(x.min(), x.max(), num=100)
            x_hat.shape = -1, 1
            y_hat = model.predict(x_hat)
            s = model.score(x, y)
            r2, corr_coef = xss(y, model.predict(x))
            # print 'R2和相关系数：', r2, corr_coef
            # print 'R2：', s, '\n'
            z = N - 1 if (d == 2) else 0
            label = '%d阶，$R^2$=%.3f' % (d, s)
            if hasattr(lin, 'l1_ratio_'):
                label += '，L1 ratio=%.2f' % lin.l1_ratio_
            plt.plot(x_hat, y_hat, color=clrs[i], lw=line_width[i], alpha=0.75, label=label, zorder=z)
        plt.legend(loc='upper left')
        plt.grid(True)
        plt.title(titles[t], fontsize=18)
        plt.xlabel('X', fontsize=16)
        plt.ylabel('Y', fontsize=16)
    plt.tight_layout(1, rect=(0, 0, 1, 0.95))
    plt.suptitle('多项式曲线拟合比较', fontsize=22)
    plt.show()

    y_max = max(max(tss_list), max(ess_rss_list)) * 1.05
    plt.figure(figsize=(9, 7), facecolor='w')
    t = np.arange(len(tss_list))
    plt.plot(t, tss_list, 'ro-', lw=2, label='TSS(Total Sum of Squares)', mec='k')
    plt.plot(t, ess_list, 'mo-', lw=1, label='ESS(Explained Sum of Squares)', mec='k')
    plt.plot(t, rss_list, 'bo-', lw=1, label='RSS(Residual Sum of Squares)', mec='k')
    plt.plot(t, ess_rss_list, 'go-', lw=2, label='ESS+RSS', mec='k')
    plt.ylim((0, y_max))
    plt.legend(loc='center right')
    plt.xlabel('实验：线性回归/Ridge/LASSO/Elastic Net', fontsize=15)
    plt.ylabel('XSS值', fontsize=15)
    plt.title('总平方和TSS=？', fontsize=18)
    plt.grid(True)
    plt.show()

结果为

线性回归：1阶，系数为： [-12.12113792   3.05477422]
线性回归：2阶，系数为： [-3.23812184 -3.36390661  0.90493645]
线性回归：3阶，系数为： [-3.90207326 -2.61163034  0.66422328  0.02290431]
线性回归：4阶，系数为： [-8.20599769  4.20778207 -2.85304163  0.73902338 -0.05008557]
线性回归：5阶，系数为： [ 21.59733285 -54.12232017  38.43116219 -12.68651476   1.98134176  -0.11572371]
线性回归：6阶，系数为： [ 14.73304785 -37.87317494  23.67462342  -6.07037979   0.42536833   0.06803132  -0.00859246]
线性回归：7阶，系数为： [ 314.30344622 -827.89446924  857.33293186 -465.46543638  144.21883851  -25.67294678    2.44658612   -0.09675941]
线性回归：8阶，系数为： [-1189.50149198  3643.69109456 -4647.92941149  3217.22814712 -1325.87384337   334.32869072   -50.57119119     4.21251817    -0.148521  ]

（2）为什么防止过拟合？

例如上述8阶的情况，得到的a0，a1等参数都非常大，那么意味着只要x发生一点点的变动，整个模型的变动就相对较大，模型很容易不稳定。

因此希望这些参数在一个合理的范围内，不需要太大也不需要太大，这样模型就会相对稳定。其中是为了平衡惩罚因子与第一项之间的比例。 

根据

的求导就能得到该公式

。

（3）其他的惩罚项来防止过拟合，L2-norm也叫作Ridge，L1-norm叫做LASSO

Ridge回归：1阶，alpha=0.166810，系数为： [-10.79755177   2.75712205]
Ridge回归：2阶，alpha=0.166810，系数为： [-2.86616277 -3.50791358  0.918189  ]
Ridge回归：3阶，alpha=0.046416，系数为： [-3.54779374 -2.8374223   0.70197829  0.02141529]
Ridge回归：4阶，alpha=0.166810，系数为： [-3.04995117 -2.03455252 -0.27692755  0.29868134 -0.02333296]
Ridge回归：5阶，alpha=0.599484，系数为： [-2.11991122 -1.79172368 -0.80325245  0.40787788 -0.02053857 -0.00138782]
Ridge回归：6阶，alpha=0.001000，系数为： [ 0.53724068 -6.00552086 -3.75961826  5.64559118 -2.21569695  0.36872911 -0.0222134 ]
Ridge回归：7阶，alpha=0.046416，系数为： [-2.3505499  -2.24317832 -1.48682673 -0.01259016  1.13568801 -0.53064855  0.09225692 -0.00560918]
Ridge回归：8阶，alpha=0.166810，系数为： [-2.11997376 -1.87289408 -1.21789779 -0.2391473   0.46199934  0.10072065 -0.11365438  0.02363517 -0.00153594]
LASSO：1阶，alpha=0.166810，系数为： [-10.83995051   2.75973604]
LASSO：2阶，alpha=0.001000，系数为： [-3.29932625 -3.31989869  0.89878903]
LASSO：3阶，alpha=0.012915，系数为： [-4.83746926 -1.48771029  0.29798381  0.05794476]
LASSO：4阶，alpha=0.012915，系数为： [-4.96121066 -1.44703055  0.32122147  0.0500902   0.00053646]
LASSO：5阶，alpha=0.166810，系数为： [-2.39014844 -2.39325347 -0.          0.16401164  0.00168607 -0.00150026]
LASSO：6阶，alpha=0.599484，系数为： [-0.         -2.9673305  -0.13745037  0.12996063  0.01895025 -0.00051774 -0.00036188]
LASSO：7阶，alpha=0.001000，系数为： [-4.51456835 -1.58477275  0.23483228  0.04900369  0.00593868  0.00044879 -0.00002625 -0.00002132]
LASSO：8阶，alpha=0.001000，系数为： [-4.62623251 -1.37717809  0.17183854  0.04307765  0.00629505  0.00069171  0.0000355  -0.00000875 -0.00000386]
ElasticNet：1阶，alpha=0.046416，l1_ratio=0.700000，系数为： [-10.86769313   2.77155741]
ElasticNet：2阶，alpha=0.046416，l1_ratio=0.100000，系数为： [-2.66468635 -3.51221793  0.9109685 ]
ElasticNet：3阶，alpha=0.012915，l1_ratio=1.000000，系数为： [-4.83746926 -1.48771029  0.29798381  0.05794476]
ElasticNet：4阶，alpha=0.012915，l1_ratio=0.100000，系数为： [-3.57655268 -2.58953187  0.55862866  0.04770441 -0.00149886]
ElasticNet：5阶，alpha=0.166810，l1_ratio=0.100000，系数为： [-1.7845673  -2.19137432 -0.26557895  0.20945337  0.00275379 -0.00193789]
ElasticNet：6阶，alpha=0.599484，l1_ratio=0.100000，系数为： [-0.84870355 -1.13574254 -1.05441649  0.26348041  0.01989297 -0.00091325 -0.00045226]
ElasticNet：7阶，alpha=0.001000，l1_ratio=0.100000，系数为： [-4.35892662 -1.71544403  0.26307018  0.04841322  0.00577139  0.00043454 -0.00002637 -0.00002108]
ElasticNet：8阶，alpha=0.001000，l1_ratio=0.950000，系数为： [-4.61710097 -1.38492703  0.17353885  0.0430437   0.00628501  0.00069073  0.00003546 -0.00000874 -0.00000386]

从上面参数的结果来看，他的确被约束了。

（4）LASSO参数降维（这里的用法我还不太懂，等后面有了新的理解来补充这部分内容）

LASSO：8阶，alpha=0.001000，系数为： [-4.62623251 -1.37717809  0.17183854  0.04307765  0.00629505  0.00069171  0.0000355  -0.00000875 -0.00000386]

上面LASSO8阶的参数中，橘色数字部分都非常的小，说明他们起到的作用也很小，那么可能不需要这么高阶的模型，此时可以将模型的参数降维。

（5）那么lamda的值该怎么取呢？

因为取不同的lamda得到的theta是不同的，那么就需要一个确定lamda的方式。

lamda不能直接算出来，但是可以选择出来：首先规定一系列的lamda值（例如在0.001-100之间取10000个数），假设使用Ridge的方式作为惩罚项，那么根据

，依次取一个lamda，同时将训练集的x和y带入，得到theta，然后根据这个theta，带入验证集中x，求得预测值y'，计算真实值y与y'的均方误差和。此时每一个lamda都能对应一个均方误差和，选取最小均方误差和对应的lamda即可。 

1.2.4 梯度下降法求解theta（第二种方式）

1.2.4.1 梯度方向 

1.2.4.2 梯度下降算法的三种参数更新方式

（1）批量梯度下降法

对所有样本计算梯度，然后对theta进行一次更新。

（2）随机梯度下降法

一个样本就算计算一次梯度，然后对theta进行一次更新。

  

（3）mini-batch梯度下降算法（这是我们最常用的方式，也是可操作性最强的方式）

选择小批量的样本计算梯度，然后对theta进行一次更新。

2. 二分类：logistic回归

那么线性模型能做分类问题吗？

答案：能，但是效果有限，如下图所示。

 

2.1 Logistic回归公式与求导 

 

2.2 Logistic回归参数估计

由于关注的问题中y值只能取0和1，自然就能想到0-1分布（这里是概率分布函数，上面的高斯函数用的是概率密度函数有所区别如何理解概率分布函数和概率密度函数？_mjiansun的博客-CSDN博客_概率分布与概率密度）。因此就有了如下假定：

那么此时整个事件的联合概率分布函数如下：

2.3 对数似然函数

相应的导数

2.4 梯度下降法参数更新

Logistic回归参数的学习规则：

比较上面的结果和线性回归的结论的差别：

它们具有相同的形式！

 

2.5 对数线性模型

一个事件的几率odds，是指该事件发生的概率与该事件不发生的概率的比值。

对数几率：logit函数

2.6 Logistic回归的损失函数

希望似然值越大越好，那么损失是希望越小越好，那么加个负号就能解决该问题。

 

 

2.7 Logistic拟合曲线

原始假设就3个特征：

X=(1,x1,x2)

theta=(theta1,theta2,theta3)

将X变为：

X=(1,x1,x2,x1^2,x2^2,x1*x2) theta=(theta1,theta2,theta3,theta4,theta5,theta6)

这样就升成6个特征维度。此时就提升了关于x1和x2的维度，但是theta的维度没有发生变化。

下图中右上角数字表示阶数。

2.8 Logistic解决异或

如下图所示，y可以取0,1值，通过x1和x2可以使用Python实现PolynomialFeatures多项式的方法来实现升维，升维后为什么能拟合任意形状曲线，我还暂时不知道。

2.9 广义线性模型（不太懂，没怎么接触过）

3. 多分类：Softmax回归（我学习到后面再来补充他）