统计学习方法（3）逻辑回归

1、从线性回归到逻辑回归（模型）

1.1 线性模型：

给定数据集{(x1,y1),(x2,y2),...,(xn,yn)}，求参数ω满足如下回归模型：
y^=ω0+ω1x1+ω2x2+...+ωnxn
线性回归模型的向量形式：
y^=ωT⋅x
其中：θT表示向量θ的转置（行向量变为了列向量）;
x为每个样本中特征值的向量形式，包括 x1，…， xn，而且 x0 =1。

模型求解：
基于最小二乘思想（方差最小）求解模型参数ω：
ω^∗=argminm1i=1∑m(ω^Tx(i)−y(i))2
写成向量形式：
ω^∗=argmin(ω^Tx−y)T(ω^Tx−y)
另Eω^=(ω^x−y)T(ω^x−y)，对其求导：
∂ω^∂Eω^=2X(XTω−y)
令导数为0，求正态方程得到：
ω^=(XTX)−1XTy

若无法求解正态方程或数据量较大，可采用梯度下降法：
对任意一点x(i)的梯度：
∂ω^∂Eω^=2x(i)(ωTx(i)−y)
随机梯度下降：
ωj+1:=ωj−α∇f(x)=ωj+α(y(i)−ωTx(i))x(i)

1.2 逻辑回归：

考虑二分类问题，由于线性回归输出值是{−∞,+∞}，应选取合适的映射函数g−(⋅)，将模型输出映射为0/1值，这里想到单位阶跃函数可以做到，但是单位阶跃函数不连续，于是想到使用类似单位阶跃函数，又能单调可微的函数，所以可以使用对数几率函数，即sigmoid函数：
y=1+e−z1
统计学习方法（3）逻辑回归-LMLPHP

线性模型带入sigmoid函数，得到：
y=1+e−(ωTx+b)1
从而得到：
log1−yy=ωTx+b
将y类后验概率估计p(y=1∣x)，则上式重写为：
logp(y=0∣x)p(y=1∣x)=ωTx+b
显然：
p(y=1∣x)=1+e(ωTx+b)e(ωTx+b)p(y=0∣x)=1+e(ωTx+b)1
为了方便，将权重向量和输入向量加以扩充，于是得到逻辑回归模型：
p(y=1∣x)=1+eωTxeωTxp(y=0∣x)=1+eωTx1
其中，ω=(ω(1),ω(2),...,ω(n),b)T，x=(x(1),x(2),...,x(n),1)T

2、策略

得到逻辑回归模型之后，考虑使用什么策略来学习模型。

2.1 极大似然估计

对于二分类，由于模型yi∈{0,1}服从0-1分布，设：
P(Y=1∣x)=ϕ(z)，P(Y=0∣x)=1−ϕ(z)则概率分布可写为：
p(yi∣xi;ω^)=ϕ(zi)yi(1−ϕ(zi))1−yi
由于概率分布已知，根据观测数据，可采用极大似然估计法来估计参数ω和b。

通过概率分布写出联合概率密度：
L(ω)=i=1∏np(y(i)∣θ,ω)=i=1∏nϕ(z(i))y(i)(1−ϕ(z(i)))(1−y(i))取对数可得：
l(ω)=i=1∑nlnp(y(i)∣θ,ω)=i=1∑ny(i)ln[ϕ(z(i))]+(1−y(i))ln[(1−ϕ(z(i)))]
对L(ω)求极大值，可求出最有可能的ω。

2.2 损失函数

同其他统计学习方法一样，逻辑回归也可以使用损失函数来学习模型，逻辑回归使用的损失函数为对数损失函数：
L(Y,P(Y∣X))=−logP(Y∣X)
该模型的损失函数为：
P(ϕ(z),y∣x;ω)=⎩⎨⎧−lnϕ(z),−ln(1−ϕ(z)),if y=1if y=0
由损失函数得到经验风险函数：
J(ω^)=−n1i=1∑n(yilnϕ(z(i))+(1−yi)ln(1−ϕ(z(i))))

3、算法

由于对数似然函数是凸函数，故可采用数值优化算法如梯度下降法、拟牛顿法求其最优解。

以下采用梯度下降法：
令
J(θ)=−l(ω)=−i=1∑ny(i)ln[ϕ(z(i))]+(1−y(i))ln[(1−ϕ(z(i)))]
∂θj∂J(θ)=−i=1∑n[y(i)ϕ(z(i))1−(1−y(i))(1−ϕ(z(i)))1]∂θj∂ϕ(z(i))
由于对于sigmoid函数：
ϕ′(z)=ϕ(z)(1−ϕ(z))
∂ω∂(ωTx+b)=x
故：
=−i=1∑n[y(i)ϕ(z(i))1−(1−y(i))(1−ϕ(z(i)))1]ϕ(z(i))(1−ϕ(z(i)))∂θj∂z(i)
=−i=1∑n[y(i)(1−ϕ(z(i)))−(1−y(i))ϕ(z(i))]xj(i)
从而得到批量梯度下降：
θj:=θj+ηi=1∑n(y(i)−ϕ(z(i)))xj(i)
随机梯度下降：
θj:=θj+η(y(i)−ϕ(z(i)))xj(i)，foriinrange(n)

4、实现

批量梯度下降：

    def train(self, X, y):
        """
        批量梯度下降：训练模型

        :param X: 输入集，n*m(n行m列)
        :param y: 标签集，1*m(1行m列)
        :return: 训练好的模型
        """
        n, m = np.shape(X)
        # X = (X, 1), w=(w;b)
        # w^T*x+b => w^T*X
        X = np.mat(np.concatenate((X, np.ones((n, 1))), axis=1))
        y = np.mat(y).T
        w = np.ones((m + 1, 1))

        alpha = self.alpha
        tol = self.tol
        records = []

        for i in range(self.max_iter):
            # 实际值与预测值之差
            error = y - 1.0 / (1 + np.exp(-X * w))
            # 误差变化记录
            records.append((error.T * error)[0, 0])
            if records[-1] <= tol:
                break
            # 负梯度方向更新参数
            w += alpha * X.T * error

        self.w = w[:-1]
        self.b = w[-1]

随机梯度下降：

    def train(self, X, y):
        """
        随机梯度下降

        :param X: 训练集
        :param y: 标签
        :return: w,b
        """
        n, m = np.shape(X)
        w = np.ones((m + 1, 1))

        alpha = self.alpha
        records = []

        for i in range(self.max_iter):
            row = random.randint(0, n - 1)
            x_i = np.mat(np.concatenate((X[row], np.ones(1))))
            y_i = y[row]
            error = y_i - 1.0 / (1 + np.exp(-np.matmul(x_i, w))[0, 0])
            records.append(error * error)

            w += alpha * error * x_i.T

        self.w = w[:-1]
        self.b = w[-1]

小蜜蜂