LogisticRegression/LinearRegression 原理及算法

线性回归

线性回归的主要思想就是通过历史数据拟合出一条直线，用这条直线对新的数据进行预测，其公式如下：
$$y^{(i)}={\theta }^T{x}^{(i)}+{\epsilon }^{(i)}$$

这里的 ε 也就是模型和实际数值之间的误差值，根据中心极限定理(许多独立随机变量组合会符合高斯分布)，我们可以接着假设误差项符合高斯分布：
$${\epsilon }^{(i)}\sim N(0,{\sigma }^2)$$

即概率密度函数为 :
$$P({\epsilon }^{(i)})=\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{({\epsilon }^{(i)}{)}^2}{2{\sigma }^2})$$

上述误差函数的概率密度函数服从高斯分布，可得：
$$P(y^{(i)}|x^{(i)};\theta )=\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2}{2{\sigma }^2})$$

即：$$y^{(i)}|x^{(i)};\theta \sim N({\theta }^T{x}^{(i)},{\sigma }^2)$$这里需要特别注意，我们不把θ认为是随机变量，而是有着我们未知的确定值，也就是把它看成我们需要去估计得到的值，也就是说上面的概率$P(y^{(i)}|x^{(i)};\theta )$ 意思是以θ为参数时，给定$x^{(i)}$条件下$y^{(i)}$的条件概率分布。

假设不同输入x(i)对应误差项ε(i)是独立同分布(IID:Independently and Identically Distributed：意思是条件独立，但是都服从同一均值方差的高斯分布)，则我们的模型可以用概率模型定义为一个极大似然估计问题：

$$L(\theta )=\prod _i^{m}P(y^{(i)}|x^{(i)};\theta )=\prod _i^m\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2}{2{\sigma }^2})$$

为方便求导变为对数似然：$$l(\theta )=lnL(\theta )=ln\prod _i^m\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2}{2{\sigma }^2})$$
$$=mln\frac{1}{\sqrt{2\pi }\sigma }+\sum _i^m-\frac{(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2}{2{\sigma }^2}$$

最大似然估计就是要求得使 $l(\theta )$ 取最大值时的 $\theta$，即使代价函数$J(\theta )$最小时的 $\theta$：$$J(\theta )=\frac{1}{2m}\sum _i^m(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2（最小二乘法）$$

Python实现

参考《机器学习算法 朱塞佩》

import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import minimize


# For reproducibility
np.random.seed(1000)

# Number of samples
nb_samples = 200


def loss(v):
    e = 0.0
    for i in range(nb_samples):
        e += np.square(v[0] + v[1]*X[i] - Y[i])
    return 0.5 * e


def gradient(v):
    g = np.zeros(shape=2)
    for i in range(nb_samples):
        g[0] += (v[0] + v[1]*X[i] - Y[i])
        g[1] += ((v[0] + v[1]*X[i] - Y[i]) * X[i])
    return g


def show_dataset(X, Y):
    fig, ax = plt.subplots(1, 1, figsize=(10, 6))

    ax.scatter(X, Y)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.grid()

    plt.show()


if __name__ == '__main__':
    # Create dataset
    X = np.arange(-5, 5, 0.05)

    Y = X + 2
    Y += np.random.uniform(-0.5, 0.5, size=nb_samples)

    # Show the dataset
    show_dataset(X, Y)

    # Minimize loss function
    result = minimize(fun=loss, x0=np.array([0.0, 0.0]), jac=gradient, method='L-BFGS-B')

    print('Interpolating rect:')
    print('y = %.2fx + %2.f' % (result.x[1], result.x[0]))

    # Compute the absolute error
    err = 0.0

    for i in range(nb_samples):
        err += np.abs(Y[i] - (result.x[1]*X[i] + result.x[0]))

    print('Absolute error: %.2f' % err)

运行结果:
Interpolating rect:
y = 1.00x +  2
Absolute error: 51.43

逻辑回归

为什么用于线性的代价函数不能用于逻辑回归?

线性回归以均方误差为代价函数。如果将其用于逻辑回归，则为参数的非凸函数。只有当函数为凸函数时，梯度下降才收敛到全局最小值。

LogisticRegression 代价函数（交叉熵函数）

Logistic Regression其思想也是基于线性回归（Logistic Regression属于广义线性回归模型）。

线性回归的公式如下：
$$z={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+{\theta }_3{x}_3...+{\theta }_n{x}_n={\theta }^{T}x$$Sigmoid 的输出为：$$h_{\theta }(x)=g(z)=\frac{1}{1+e^{-z}}=\frac{1}{1+e^{-{\theta }^{T}x}}$$是介于（0，1）之间的，中间值是0.5，这就表明了数据属于某一类别的概率，例如 ：
$h_{\theta }(x)$ < 0.5 则说明当前数据属于A类；
$h_{\theta }(x)$ > 0.5 则说明当前数据属于B类。
所以我们可以将sigmoid函数看成样本数据的概率密度函数。

交叉熵损失函数:

交叉熵损失函数的推导 ：

同线性回归，我们可以把模型用概率表示：

$$P(y^{(i)}=1|x^{(i)};\theta )=h_{\theta }(x^{(i)})$$

$$P(y^{(i)}=0|x^{(i)};\theta )=1-h_{\theta }(x^{(i)})$$

我们可以进一步把两式整合：

$$P(y^{(i)}|x^{(i)};\theta )=h_{\theta }(x^{(i)}{)}^{y^{(i)}}(1-h_{\theta }(x^{(i)}){)}^{(1-y^{(i)})}$$

同样我们可以把模型最优问题看做是极大似然估计问题：

$$L(\theta )=\prod _i^{m}P(y^{(i)}|x^{(i)};\theta )=\prod _i^m{h}_{\theta }(x^{(i)}{)}^{y^{(i)}}(1-h_{\theta }(x^{(i)}){)}^{(1-y^{(i)})}$$

取对数似然：

$$lnL(\theta )=\sum _i^m\left(y^{(i)}ln(h_{\theta }(x^{(i)}))+(1-y^{(i)})ln(1-h_{\theta }(x^{(i)}))\right)$$

最大似然估计就是要求得使$lnL(\theta )$取最大值时的 $\theta$，即使代价函数$J(\theta )$最小时的 $\theta$ ，这里可以使用梯度上升法求解。我们稍微变换一下：

$$J(\theta )=-\frac{1}{m}\sum _i^m\left(y^{(i)}ln(h_{\theta }(x^{(i)}))+(1-y^{(i)})ln(1-h_{\theta }(x^{(i)}))\right)$$

因为乘了一个负的系数$-\frac{1}{m}$，然后就可以使用梯度下降算法进行参数求解了。

梯度下降法求J(θ)的最小值

求J(θ)的最小值可以使用梯度下降法，根据梯度下降法可得θ的更新过程

$${\theta }_j={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}J(\theta ),(j=0,1,...,n)$$

以下步骤参考：Logistic回归总结 作者：洞庭之子

Python实现

def weightInitialization(n_features):
    w = np.zeros((1,n_features))
    b = 0
    return w,b
    
def sigmoid_activation(result):
    final_result = 1/(1+np.exp(-result))
    return final_result

def model_optimize(w, b, X, Y):
    m = X.shape[0]
    #Prediction
    final_result = sigmoid_activation(np.dot(w,X.T)+b)
    Y_T = Y.T
    cost = (-1/m)*(np.sum((Y_T*np.log(final_result)) + ((1-Y_T)*(np.log(1-final_result)))))
    #Gradient calculation
    dw = (1/m)*(np.dot(X.T, (final_result-Y.T).T))
    db = (1/m)*(np.sum(final_result-Y.T))
    grads = {"dw": dw, "db": db}
    return grads, cost
    
def model_predict(w, b, X, Y, learning_rate, no_iterations):
    costs = []
    for i in range(no_iterations):
        grads, cost = model_optimize(w,b,X,Y)
        dw = grads["dw"]
        db = grads["db"]
        #weight update
        w = w - (learning_rate * (dw.T))
        b = b - (learning_rate * db)
        if (i % 100 == 0):
            costs.append(cost)
            #print("Cost after %i iteration is %f" %(i, cost))
    #final parameters
    coeff = {"w": w, "b": b}
    gradient = {"dw": dw, "db": db}
    return coeff, gradient, costs
    
def predict(final_pred, m):
    y_pred = np.zeros((1,m))
    for i in range(final_pred.shape[1]):
        if final_pred[0][i] > 0.5:
            y_pred[0][i] = 1
    return y_pred

参考：

• [ML] 逻辑回归——详细概述 数据科学与人工智能

• 机器学习算法之：指数族分布与广义线性模型

为了预测数据点属于哪个类别，可以设置一个阈值。根据这个阈值，将获得的估计概率划分为类别。如果predicted_value ≥ 0.5，电子邮件邮件分类为垃圾邮件，反之不是。

决策边界可以是线性的，也可以是非线性的。多项式增加阶数可获得复杂的决策边界。

逻辑回归的模型本质上是一个线性回归模型，逻辑回归都是以线性回归为理论支持的。但线性回归模型无法做到sigmoid的非线性形式，sigmoid可以轻松处理0/1分类问题。