逻辑回归学习及实现

##逻辑回归原理

将样本的特征与样本的概率联系起来，通过概率来判定所属分类；先构建多元线性方程y，然后将y带入阶跃函数(step function)sigmoid，转化成求分类概率的问题(一般概率大于0.5分类为1，小于0.5分类为0),而Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化 算法来完成。在最优化算法中，最常用的就是梯度上升(下降)算法，而梯度上升(下降)算法又可以简化为随机 梯度上升(下降)算法。。

注：通常逻辑回归会用作二分类，也可以实现多分类。

##逻辑回归的优缺点

优点：

1. 预测结果是界于0和1之间的概率；
2. 可以适用于连续性和类别性自变量；
3. 计算代价不高，易于理解和实现。

缺点：

1. 对模型中自变量多重共线性较为敏感，例如两个高度相关自变量同时放入模型，可能导致较弱的一个自变量回归符号不符合预期，符号被扭转。​需要利用因子分析或者变量聚类分析等手段来选择代表性的自变量，以减少候选变量之间的相关性。（解决多重共线性的几种方法：https://blog.csdn.net/nieson2012/article/details/48980491；https://blog.csdn.net/diyiziran/article/details/17025471）
2. 预测结果呈“S”型，因此从log(odds)向概率转化的过程是非线性的，在两端随着​log(odds)值的变化，概率变化很小，边际值太小，slope太小，而中间概率的变化很大，很敏感。 导致很多区间的变量变化对目标概率的影响没有区分度，无法确定阀值。(https://blog.csdn.net/yezi_1026/article/details/53121723)

逻辑回归实现步骤

1. 寻找h函数（即hypothesis）；

2. 构造J函数（损失函数）；
3. 使得J函数最小并求得回归参数（θ）

构造预测函数$h_{\theta }(x)$：

$$g(z)=\frac{1}{1+e^{-z}}$$
其中$x_0=1$, $z$的表达式：
$$z={\theta }_0{x}_0+{\theta }_1{X}_1+{\theta }_2{x}_2+...+{\theta }_n{x}_n=\sum _{i=0}^n{\theta }_i{x}_i={\theta }^T\bullet x$$
$$h_{\theta }(x)=g({\theta }^T\bullet x)=\frac{1}{1+e^{-{\theta }^T\bullet x}}$$

构造损失函数$J(\theta )$

这里将$h_{\theta}(x) $对数化：

这里我们可以对$J(\theta )$构建似然函数。

构建似然函数：函数$h_{\theta }(x)$值有特殊含义，表示取1的概率，因此对于输入x分类结果为类别1和类别0的概率分别为
$$\begin{gathered} P(y=1|x,\theta )=h_{\theta }(x) \\ P(y=0|x,\theta )=1-h_{\theta }(x) \end{gathered}$$
进一步可以推导成：
$$P(y=1|x,\theta )=(h_{\theta }(x){)}^y(1-h_{\theta }(x){)}^{1-y}$$
取似然函数：
$$L(\theta )=\prod _{i=1}^m(h_{\theta }(x){)}^y(1-h_{\theta }(x){)}^{1-y}$$
对$L(\theta )$两边取对数：

最大似然估计就是求使取最大值时的θ，其实这里可以使用梯度上升法求解，求得的θ就是要求的最佳参数。但是，在Andrew Ng的课程中将取为下式，即：
$$J(\theta )=-\frac{1}{m}l(\theta )$$
因为乘了一个负的系数-1/m，所以取最小值时的θ为要求的最佳参数。

使用梯度下降法逼近$J(\theta )$的最小值：

$\nabla J(\theta )$的更新式子：
$$\nabla J(\theta )=\frac{\partial J(\theta )}{\partial {\theta }_j}$$

整理得到：$$\nabla J(\theta )=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^i)-y^i)x_j^i$$
最后向量化$\nabla J(\theta )$得到：
$$\nabla J(\theta )=\frac{1}{m}{X}_b^T\bullet (g(X_b\bullet \theta )-y)$$
$注：j表示维度，i表示第i个样本。$

代码封装实现

import numpy as np
from sklearn.metrics import accuracy_score

class LogisticRegression:

    def __init__(self):
        """初始化Logistic Regression模型"""
        self.coef_ = None
        self.intercept_ = None
        self._theta = None

    def _sigmoid(self, t):
        return 1. / (1 + np.exp(-t))
    
    def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
        """根据训练数据集X_train, y_train, 使用梯度下降法训练Logistic Regression模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"
        

        def J(theta, X_b, y):
            y_hat = self._sigmoid(X_b.dot(theta))
            
            try:
                return -np.sum(y*np.log(y_hat) + (1-y)*(1-np.log(1-y_hat))) / len(y)
            except:
                return float('inf')

        def dJ(theta, X_b, y):
            return X_b.T.dot(self._sigmoid(X_b.dot(theta)) - y)  / len(X_b)

        def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):

            theta = initial_theta
            cur_iter = 0

            while cur_iter < n_iters:
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient
                if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
                    break

                cur_iter += 1

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters)

        self.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    
    def predict_proba(self, X_predict):
        """给定待预测数据集X_predict，返回表示X_predict的概率向量"""
        assert self.intercept_ is not None and self.coef_ is not None, \
            "must fit before predict!"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"

        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        return self._sigmoid(X_b.dot(self._theta))
    
    def predict(self, X_predict):
        """给定待预测数据集X_predict，返回表示X_predict的结果向量"""
        assert self.intercept_ is not None and self.coef_ is not None, \
            "must fit before predict!"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"

        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        proba = self.predict_proba(X_predict)
        return np.array(proba>=0.5, dtype='int')

    def score(self, X_test, y_test):
        """根据测试数据集 X_test 和 y_test 确定当前模型的准确度"""

        y_predict = self.predict(X_test)
        return accuracy_score(y_test, y_predict)

    def __repr__(self):
        return "LogisticRegression"

验证封装代码

from sklearn.datasets import load_iris
from matplotlib.pyplot as plt
X = load_iris().data
y = load_iris().target
X = X[y<2,:2]
y = y[y<2]
plt.scatter(X[y==0,0],X[y==0,1],color= 'r')
plt.scatter(X[y==1,0],X[y==1,1],color = 'g')

Logr = LogisticRegression()
Logr.fit_gd(X,y)
Logr.predict_proba(X)
Logr.predict(X)
Logr.score(X,y)

0.98999999999999999

学习笔记参考：
《机器学习实战》和《Python3入门机器学习 经典算法与应用》
优缺点内容参考：https://blog.csdn.net/mr_hhh/article/details/79433094
似然估计参考：https://blog.csdn.net/xiaoshunzi111/article/details/52329971