《机器学习》线性模型公式推导与算法实现

线性回归

参考西瓜书《机器学习》线性回归

给定训练集D={(\boldsymbol x_1, y_1), (\boldsymbol x_2, y_2), ..., (\boldsymbol x_i, y_i), ( \boldsymbol x_n, y_n)},其中\boldsymbol {x_i} = (x_{i1};x_{i2}; ...; x_{im})y_i\in \mathbb{R}.线性回归(linear regression)试图学得一个线性模型以尽可能准确地预测实值输出标记

以一元线性回归为例,来详细讲解wb的最小二乘法估计。

线性回归试图学得:f(x_i)=wx_i+b,使得f(x_i)\simeq y_i

  • 最小二乘法是基于预测值和真实值的均方差最小化的方法来估计参数wb,即:

    \eqalign{ (w^*,b^*) & = \mathop {\arg \min }\limits_{(w,b)} \sum\limits_{i = 1}^n {{{(f({x_i}) - {y_i})}^2}} \cr & = \mathop {\arg \min }\limits_{(w,b)} \sum\limits_{i = 1}^n {{{({y_i} - w{x_i} - b)}^2}} \cr}

  • 求解wb使{E_{(w,b)}} = \sum\limits_{i = 1}^n {{{({y_i} - w{x_i} - b)}^2}}最小化的过程,称为线性回归模型的最小二乘“参数估计”(parameter estimation)。
  • 我们可将E_{(w,b)}分别对wb求导,得到

    \eqalign{ & {{\partial E(w,b)} \over {\partial w}} = 2\left( {w\sum\limits_{i = 1}^n {x_i^2 - \sum\limits_{i = 1}^n {({y_i} - b){x_i}} } } \right) \cr & {{\partial E(w,b)} \over {\partial b}} = 2\left( {nb - \sum\limits_{i = 1}^n {({y_i} - w{x_i})} } \right) \cr}

  • 令上两式为零可得到wb最优解的闭式(closed-form)解:

    \eqalign{ w & = {{\sum\limits_{i = 1}^n {{y_i}({x_i} - \bar x)} } \over {\sum\limits_{i = 1}^n {x_i^2 - {\textstyle{1 \over n}}{{\left( {\sum\limits_{i = 1}^n {{x_i}} } \right)}^2}} }} \cr b & = {\textstyle{1 \over n}}\sum\limits_{i = 1}^n {({y_i} - w{x_i})} \cr}

    其中,\bar x = {\textstyle{1 \over n}}\sum\limits_{i = 1}^n {{x_i}}x的均值。

更一般的情况

给定数据集D,样本由d个属性描述,此时我们试图学得

f(\boldsymbol {x_i}) = \boldsymbol {w^T}\boldsymbol {x_i} + b, 使得f({x_i}) \simeq {y_i}

这称为“多元线性回归”(multivariate linear regression).

类似的,可利用最小二乘法来对wb进行估计。为便于讨论,我们把wb吸收入向量形式\widehat w = (w;b)

相应的,把数据集D表示为一个mx(d+1)大小的矩阵{\bf{X}},其中每行对应于一个示例,该行前d个元素对应于前d个属性值,最后一个元素恒置为1,即
{\bf{X}} = \left( {\matrix{ {{x_{11}}} & \cdots & {{x_{1d}}} & 1 \cr {{x_{21}}} & \cdots & {{x_{2d}}} & 1 \cr \vdots & \ddots & \vdots & 1 \cr {{x_{m1}}} & \cdots & {{x_{md}}} & 1 \cr } } \right) = \left( {\matrix{ {x_1^T} \cr {x_2^T} \cr \vdots \cr {x_m^T} \cr } \matrix{ 1 \cr 1 \cr \vdots \cr 1 \cr } } \right)

再把标记也写成向量形式\boldsymbol{y}=(y_1;y_2;...;y_m),有

\hat {\boldsymbol {w}^*} = \mathop {\arg \min }\limits_{\hat w} \boldsymbol{(y - X\hat w)^T}\boldsymbol {(y - X\hat w)}

{E_{\hat w}} = {(\boldsymbol y - {\bf{X}}\boldsymbol {\hat w})^T}(\boldsymbol y - {\bf{X}}\boldsymbol {\hat w}),对\hat w求导得到:

{{\partial {E_{\hat w}}} \over {\partial \hat {\boldsymbol w}}} = 2{{\bf{X}}^T}({\bf{X}}\hat w - \boldsymbol y)

令上式为零可得\hat {w}最优解的闭式解。

{{\bf{X}}^T}{\bf{X}}为满秩矩阵或正定矩阵时,令上式为零可得:
\boldsymbol {\hat w^*} = {({{\bf{X}}^T}{\bf{X}})^{ - 1}}{{\bf{X}}^T}\boldsymbol {y}
\boldsymbol {\hat x_i} = (\boldsymbol{x_i};1),则最终学得的多元线性回归模型为
f({\boldsymbol{\hat x_i}}) = \boldsymbol {\hat x_i}{({{\bf{X}}^T}{\bf{X}})^{ - 1}}{{\bf{X}}^T}\boldsymbol y

广义线性模型

线性回归假定输入空间到输出空间的函数映射成线性关系,但现实应用中,很多问题都是非线性的。为了拓展其应用场景,我们可以将线性回归的预测值做一个非线性的函数变化去逼近真实值,这样得到的模型统称为广义线性模型(generalized linear model):
y = {g^{ - 1}}(\boldsymbol{w^T}\boldsymbol x + b)
其中函数g(\cdot)称为“联系函数”(link function),它连续且充分光滑

当联系函数为g(\cdot)=ln(\cdot)时,称为对数线性回归。即
\ln y= \boldsymbol{w^T}\boldsymbol x + b \\ y = e^{\boldsymbol{w^T}\boldsymbol x + b}

逻辑斯蒂回归

前面介绍的都是完成回归任务,如果要做的是分类任务该怎么办呢?

按照上面介绍的广义线性模型,要完成分类任务,也就是去寻找一个合适的g(\cdot)函数。为了简化问题,我们先考虑二分类任务,其输出标记y \in \{ 0,1\}。线性模型产生的预测值z=\boldsymbol{w^T}\boldsymbol {x} + b是实值,因此,我们需要将实值z转换为0/1值。最理想的是“单位阶跃函数”。
y=\begin{cases} 0,&z<0 \\ 0.5, &z=0 \\ 1,&z>0 \end{cases}
即若预测值z大于零就判为正例,小于零则判为反例,预测值为临界值零则可任意判定,如图所示

单位阶跃函数与对数几率函数

从图中可以发现,单位阶跃函数不连续,因此不能直接用作联系函数g(\cdot)。于是我们希望找到能在一定程度上近似单位阶跃函数的替代函数,并希望它在临界点连续且单调可微。

对数几率函数(logistic function)正是这样一个常用的替代函数。
y = \frac{1}{{1 + {e^{ - (\boldsymbol{w^T}\boldsymbol x + b)}}}}
从图中可以看出,它的图像形式S,对这类图像形式S的函数称为"Sigmoid函数",对数几率函数是它最重要的代表。

这种利用对数几率函数去逼近真实标记的对数几率的模型称为“对数几率回归”,又常叫做“逻辑斯蒂回归”,虽然它的名字是“回归”,但实际上是一种分类学习方法。

对数逻辑回归有很多优点:

  • 1.可以直接对分类可能性进行预测,将y视为样本x作为正例的概率;
  • 2.无需事先假设数据分布,这样避免了假设分布不准确带来的问题;
  • 3.它是任意阶可导的凸函数,可直接应用现有的数值优化算法求取最优解。

参数\boldsymbol wb的确定

将y视为样本\boldsymbol x作为正例的可能性。显然有:
p(y=1|\boldsymbol x) = \frac{{{e^{\boldsymbol {w^T}\boldsymbol x + b}}}}{{1 + {e^{\boldsymbol {w^T}\boldsymbol x + b}}}}
p(y=0|\boldsymbol x) = \frac{1}{{1 + {e^{\boldsymbol {w^T}\boldsymbol x + b}}}}

于是,可通过“极大似然法”来估计参数\boldsymbol wb。给定数据集\{ (\boldsymbol {x_i},{y_i})\} _{i = 1}^m,其“对数似然”为:
l(\boldsymbol w,b) = \sum\limits_{i = 1}^n {\ln p({y_i}|\boldsymbol {x_i};\boldsymbol w,b)}

\eqalign{ & \boldsymbol \beta = (\boldsymbol w;b),\hat {\boldsymbol x} = (\boldsymbol x;1) \cr & {\boldsymbol \beta ^T}\hat {\boldsymbol x} = \boldsymbol{w^T}\boldsymbol x + b \cr & {p_1}(\boldsymbol {\hat x};\boldsymbol \beta ) = p(y = 1|\boldsymbol {\hat x};\boldsymbol \beta ) \cr & {p_0}(\boldsymbol {\hat x};\boldsymbol \beta ) = p(y = 0|\boldsymbol {\hat x};\boldsymbol \beta ) \cr}
则似然项可重写为:
p({y_i}|\boldsymbol {x_i};\boldsymbol w,b) = {y_i}{p_1}({\boldsymbol {\hat x_i}};\boldsymbol \beta ) + (1 - {y_i}){p_0}(\boldsymbol {\hat x};\boldsymbol \beta )
将上式带入似然函数:

\eqalign{ l(\beta ) & = \sum\limits_{i = 1}^n {\ln ({y_i}{p_1}(\hat x;\beta ) + (1 - {y_i}){p_0}(\hat x;\beta )} \cr & = \sum\limits_{i = 1}^n {\ln \left( {{y_i}\frac{{{e^{{\beta ^T}\hat x}}}}{{1 + {e^{{\beta ^T}\hat x}}}} + (1 - {y_i})\frac{1}{{1 + {e^{{\beta ^T}\hat x}}}}} \right)} \cr & = \sum\limits_{i = 1}^n {\ln \left( {\frac{{{y_i}{e^{{\beta ^T}\hat x}} + (1 - {y_i})}}{{1 + {e^{{\beta ^T}\hat x}}}}} \right)} \cr & = \sum\limits_{i = 1}^n {\ln \left( {\frac{{{y_i}({e^{{\beta ^T}\hat x}} - 1) + 1}}{{1 + {e^{{\beta ^T}\hat x}}}}} \right)} \cr & = \sum\limits_{i = 1}^n {\left( {\ln ({y_i}({e^{{\beta ^T}\hat x}} - 1) + 1) - \ln (1 + {e^{{\beta ^T}\hat x}})} \right)} \cr & =\begin{cases} {\sum\limits_{i = 1}^n { - \ln (1 + {e^{{\beta ^T}\hat x}})} },&y_i=0 \\ {\sum\limits_{i = 1}^n {({\beta ^T}\hat x - \ln (1 + {e^{{\beta ^T}\hat x}}))} },&y_i=1 \cr \end{cases}}
考虑到y_i \in \{0, 1\},即最大化l(\boldsymbol w,b)等价于最小化
l(\boldsymbol \beta ) = \sum\limits_{i = 1}^n {(-y_i\boldsymbol{\beta ^T}\boldsymbol {\hat x_i} + \ln (1 + {e^{{\beta ^T}\hat x_i}}))}

接下来可以利用数值优化算法对其求解,即
\beta ^* {\text{ = }}\mathop {\arg \min }\limits_\beta l(\beta )

对数逻辑回归代码实现

回归模型
y = \frac{1}{{1 + {e^{ - z}}}}
损失函数(最小化):
l(\boldsymbol \beta ) = \sum\limits_{i = 1}^n {(-y_i\boldsymbol{\beta ^T}\boldsymbol {\hat x_i} + \ln (1 + {e^{{\beta ^T}\hat x_i}}))}

\beta ^* {\text{ = }}\mathop {\arg \min }\limits_\beta l(\beta )
以牛顿法求解为例:其第t+1轮迭代解的更新公式为
\boldsymbol \beta ^{t+1}=\boldsymbol \beta ^t-\left ( \frac{\partial ^2l(\boldsymbol \beta))}{\partial {\boldsymbol \beta}\partial {\boldsymbol \beta} ^T} \right )^{-1}\frac{\partial l(\boldsymbol \beta)}{\partial \boldsymbol \beta}
其中关于\boldsymbol \beta的一阶、二阶导数分别为
\frac{\partial l(\beta)}{\partial \beta} = -\sum ^m_{i=1}\hat x_i(y_i-p_1(\hat x_i;\beta))
\frac{\partial ^2l(\beta)}{\partial \beta\partial \beta ^T}=\sum ^m_{i=1}\hat x_i\hat x_i^Tp_1(\hat x_i;\beta)(1-p_1(\hat x_i; \beta))

import os
import numpy as np
import pandas as pd
data = pd.read_csv("../data/irisdata.txt")
# 只保留两种标签,进行二分类任务
data = data[data['name'] != 'Iris-setosa']
data['name'].value_counts()
Iris-versicolor    50
Iris-virginica     50
Name: name, dtype: int64
# 分离标签,并将标签映射到数值
y = data['name']
y[y == 'Iris-versicolor'] = 1
y[y == 'Iris-virginica'] = 0
X = data.drop('name', axis=1)
# 划分训练集和验证集
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5)
class LogisticReressionClassifier:
    def __init__(self, max_iter):
        self.max_iter = max_iter
        
    def sigmod(self, z):
        return 1 / (1 + np.exp(-z))
    
    def fit(self, X, y):
        self.beta = np.random.normal(size=(X.shape[0], X.shape[1] + 1))  # 初始化参数
        self.X_hat = np.c_[X, np.ones(X.shape[0])]  # 为数据集加入一列1
        self.loss_function(X, y)  # 打印训练前loss
        for j in range(self.max_iter): # 迭代优化
            pd1 = 0  # 一阶偏导
            for i in  range(len(y)):
                pd1 -= self.X_hat[i]*(y[i] - self.sigmod(np.dot(self.beta[i].T, self.X_hat[i])))
            pd2 = 0  # 二阶偏导
            for i in range(len(y)):
                pd2 += self.X_hat[i].dot(self.X_hat[i].T.dot(self.sigmod(self.beta[i].T.dot(self.X_hat[i]))*(1 - self.sigmod(self.beta[i].T.dot(self.X_hat[i])))))
            self.beta = self.beta - (1 / pd2)*pd1  # 更新参数beta
        self.loss_function(X, y)  # 打印训练后的loss
        
    def loss_function(self, X, y):
        loss = 0
        # 根据损失函数公式计算当前loss
        for i in range(len(y)):
            loss += -y[i]*np.dot(self.beta[i].T, self.X_hat[i]) + np.log(1 + np.exp(np.dot(self.beta[i].T, self.X_hat[i])))
        print(loss)
        
    def predict(self, X):
        y = [] # 存储预测结果
        X = np.c_[X, np.ones(X.shape[0])]  # 为训练集加入一列1
        for i in range(X.shape[0]):
            # 计算样本作为正例的相对可能性(几率)
            odds = self.sigmod(np.mean(self.beta, axis=0).T.dot(X[i]))
            if (odds >= 0.5):
                y.append(1)
            else:
                y.append(0)
        return y

clf = LogisticReressionClassifier(10000)
clf.fit(X_train.values, y_train.values)
187.27577618364464
38.2785420108109
y_pred = clf.predict(X_test.values)
# 正确率
sum(y_pred == y_test)/len(y_test)
0.96

更多

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 219,110评论 6 508
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,443评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 165,474评论 0 356
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,881评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,902评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,698评论 1 305
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,418评论 3 419
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,332评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,796评论 1 316
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,968评论 3 337
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,110评论 1 351
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,792评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,455评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,003评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,130评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,348评论 3 373
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,047评论 2 355

推荐阅读更多精彩内容