[机器学习-从入门到入土] 逻辑回归

个人导航

知乎：https://www.zhihu.com/people/byzh_rc

CSDN：https://blog.csdn.net/qq_54636039

logistic-regression.ipynb

逻辑回归的假设模型为：
$$h_{\theta }(x)=g({\theta }^{\mathrm{T}}x)$$
​

其中函数 $g(\cdot )$ 是Sigmoid函数，定义为：
$$g(z)=\frac{1}{1+\exp (-z)}$$

g = 1 / (1 + np.exp(-z))

1.不带正则项

逻辑回归的代价函数为: (不带正则项)
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m[-y^{(i)}\log (h_{\theta }(x^{(i)}))-(1-y^{(i)})\log (1-h_{\theta }(x^{(i)})\left)\right]$$

h = sigmoid(X @ theta)
h = np.clip(h, 1e-10, 1 - 1e-10)
J = np.sum(
    -y * np.log(h) - (1 - y) * np.log(1 - h)
) / m
# X: (m, input_size) y: (m,)
# theta: (input_size,)
# h: (m,)
# -> J: 数

对应的梯度向量各分量为: (不带正则项)
$$\frac{\partial J(\theta )}{\partial {\theta }_j}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)}$$

h = sigmoid(X @ theta)
grad = (X.T @ (h - y)) / m
# X: (m, input_size) y: (m,)
# theta: (input_size,)
# h: (m,)
# -> grad: (input_size,)

有: $1-g=\frac{exp(-z)}{1+exp(-z)}=exp(-z)\ast g$ -> $1-h=exp(-z)\ast h$
有: $g^{\prime }=exp(-z)\ast g^2$ -> $h^{\prime }=exp(-z)\ast h^2$
则:
$$h^{\prime }=h\ast (1-h)$$
整理得:
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}\right(-\log (h_{\theta }(x^{(i)}))+log(1-h_{\theta }(x^{(i)})))-log(1-h_{\theta }(x^{(i)}))]$$
故:
$$\frac{\partial J(\theta )}{\partial z}=\frac{1}{m}\sum _{i=1}^m[y^{(i)}(-\frac{h^{\prime }}{h}+\frac{-h^{\prime }}{1-h})-\frac{-h^{\prime }}{1-h}]=\frac{1}{m}\sum _{i=1}^m[y^{(i)}(h-1-h)+h]$$
上式是对$z$求导, 而我们需要对${\theta }_j$求导
$$\frac{\partial h}{\partial {\theta }_j}=\frac{\partial h}{\partial z}\cdot \frac{\partial z}{\partial {\theta }_j}=h^{\prime }(z)\cdot x_j^{(i)}$$
故:
$$\frac{\partial J(\theta )}{\partial {\theta }_j}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)}$$

发现, 逻辑回归的代价函数(无正则项)求导后居然和线性回归一样

使用**非线性共轭梯度算法(nonlinear conjugate gradient algorithm)**来优化函数:

import scipy.optimize as op
ret = op.fmin_cg(cost_function, # 最小化目标函数
                 theta, # 初始参数
                 fprime=cost_gradient, # 提供梯度函数 
                 # 有了梯度，优化器就知道沿哪个方向下降最快，从而加快收敛
                 args=(X, y), # 其他参数通过 args 传递
                 maxiter=400, # 最大迭代次数
                 full_output=True) # 返回完整输出信息
theta_opt, cost_min, _, _, _ = ret

cost_function为代价函数
theta 为需要优化的参数初始值
fprime=cost_gradient 给出了代价函数的梯度
args=(X, y) 给出了需要优化的函数与对应的梯度计算所需要的其他参数
maxiter=400 给出了最大迭代次数

2.带正则项

特征转换:
创建更多的特征是充分挖掘数据中的信息的一种有效手段
将数据映射为其六阶多项式的所有项
$$\text{map\_feature}(\mathbf{x})=\left[\begin{array}{c}1\\ x_1\\ x_2\\ x_1^2\\ x_1{x}_2\\ x_2^2\\ x_1^3\\ ⋮\\ x_1{x}_2^5\\ x_2^6\end{array}\right]$$

$x_1^i{x}_2^j,i+j\le k(最大阶数)$

$k=0$ -> $1$种情况
$k=1$ -> $2$种情况
…
$k=d$ -> $d+1$种情况

所以k=6有7种情况, 则k从0到6, 有1+2+3+4+5+6+7=28项
-> input_size变为28

逻辑回归的代价函数为: (正则项)
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m\left[-y^{(i)}\log \left(h_{\theta }(x^{(i)})\right)-(1-y^{(i)})\log \left(1-h_{\theta }(x^{(i)})\right)\right]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$

h = sigmoid(X @ theta)
h = np.clip(h, 1e-10, 1 - 1e-10)
J += np.sum(-y * np.log(h) - (1 - y) * np.log(1 - h)) / m
J += lmb / (2 * m) * np.sum(theta[1:] ** 2)
# X: (m, input_size) y: (m,)
# theta: (input_size,)
# h: (m,)
# -> J: 数

对应的梯度向量各分量为: (正则项)
$$\begin{array}{rl}\frac{\partial J(\theta )}{\partial {\theta }_0}& =\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_0^{(i)}\text{for }j=0\\ \frac{\partial J(\theta )}{\partial {\theta }_j}& =\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)}+\frac{\lambda }{m}{\theta }_j\text{for }j\ge 1\end{array}$$

h = sigmoid(X @ theta)
grad = (X.T @ (h - y)) / m
grad[1:] = grad[1:] + (lmb / m) * theta[1:]
# X: (m, input_size) y: (m,)
# theta: (input_size,)
# h: (m,)
# -> grad: (input_size,)

python代码

import numpy as np
import scipy.optimize as op
import matplotlib.pyplot as plt

######################### cell 8 #########################

def plot_data(X, y):
    """This function plots the data points X and y into a new figure.
    It plots the data points with red + for the positive examples,
    and blue o the negative examples. X is assumed to be a Mx2 matrix.
    
    X: shape:nx2 
    y: shape:nx1
    """

    plt.figure()
    # ====================== YOUR CODE HERE ======================
    plt.scatter(X[y==1, 0], X[y==1, 1], color='r', marker='+', s=50, label='Admitted')
    plt.scatter(X[y==0, 0], X[y==0, 1], color='b', marker='o', s=50, label='Not admitted')
    # ============================================================

    plt.xlabel("Exam 1 Score")
    plt.ylabel("Exam 2 Score")
    plt.legend(loc='upper right')

######################### cell 10 #########################

# 加载数据
data = np.loadtxt("LR_data1.txt", delimiter=",")
X, y = data[:, :2], data[:, 2]
# print(X.shape) # (100, 2)
# print(y.shape) # (100,)

# 可视化数据
# ====================== YOUR CODE HERE ================
plot_data(X, y)
# ======================================================
plt.show()

######################### cell 12 #########################

def plot_decision_boundary(theta, X, y):
    """绘制分类面。"""
    plot_data(X[:, 1:], y)

    _, d = X.shape

    if d <= 3:
        plot_x = np.array([np.min(X[:, 1])-2, np.max(X[:, 1])+2])
        plot_y = -1.0 / theta[2]*(theta[1]*plot_x + theta[0])
        plt.plot(plot_x, plot_y, 'm-', label="Decision Boundary")

        plt.xlim([30, 100])
        plt.ylim([30, 100])
    else:
        n_grid = 50
        u = np.linspace(-1, 1.5, n_grid)
        v = np.linspace(-1, 1.5, n_grid)

        z = np.zeros((n_grid, n_grid))

        for i in range(n_grid):
            for j in range(n_grid):
                uu, vv = np.array([u[i]]), np.array([v[j]])
                z[i, j] = np.dot(map_feature(uu, vv), theta)

        z = z.T

        CS = plt.contour(u, v, z, linewidths=2, levels=[0.0], colors=['m'])
        CS.collections[0].set_label('Decision boundary')

    plt.legend()

######################### cell 14 #########################

def sigmoid(z):
    """Compute sigmoid function"""

    z = np.asarray(z)
    g = np.zeros_like(z)

    # ====================== YOUR CODE HERE ======================
    g = 1 / (1 + np.exp(-z))
    # ============================================================

    return g

######################### cell 15 #########################

# 测试 sigmoid 函数
z = np.array([-10.0, -5.0, 0.0, 5.0, 10.0])
g = sigmoid(z)
print("Value of sigmoid at [-10, -5, 0, 5, 10] are:\n", g)

######################### cell 17 #########################

def cost_function(theta, X, y):
    """逻辑回归的代价函数，无正则项。"""

    J = 0.0
    # ====================== YOUR CODE HERE ======================
    m = len(y)
    h = sigmoid(X@theta)
    h = np.clip(h, 1e-10, 1 - 1e-10)
    J = np.sum(
        -y*np.log(h) - (1-y)*np.log(1-h)
    ) / m
    # ============================================================
    return J

######################### cell 18 #########################

def cost_gradient(theta, X, y):
    """逻辑回归的代价函数的梯度，无正则项。"""
    m = 1.0*len(y)
    grad = np.zeros_like(theta)
    # ====================== YOUR CODE HERE ======================
    h = sigmoid(X@theta)
    grad = (X.T @ (h - y)) / m
    # ============================================================
    return grad

######################### cell 20 #########################

def predict(theta, X):
    """Predict whether the label is 0 or 1
    using learned logistic regression parameters theta.
    input： theta：model's parameters
            X: input samples
    output：0 or 1
    """
    m, _ = X.shape
    pred = np.zeros((m, 1), dtype=bool)

    # ====================== YOUR CODE HERE ======================
    pred = sigmoid(X@theta) > 0.5
    # ============================================================

    return pred

######################### cell 22 #########################

def logistic_regression():
    """针对第一组数据建立逻辑回归模型。"""

    # 加载数据
    data = np.loadtxt("LR_data1.txt", delimiter=",")
    X, y = data[:, :2], data[:, 2]

    # 计算代价与梯度
    m, _ = X.shape
    X = np.hstack((np.ones((m, 1)), X))

    # 初始化参数
    theta_initial = np.zeros_like(X[0])

    # 计算并打印初始参数对应的代价与梯度
    cost = cost_function(theta_initial, X, y)
    grad = cost_gradient(theta_initial, X, y)
    print("Cost at initial theta (zeros): ", cost)
    print("Gradient at initial theta (zeros): \n", grad)

    # 使用 scipy.optimize.fmin_cg 优化模型参数
    args = (X, y)
    maxiter = 200
    # ====================== YOUR CODE HERE ======================
    ret = op.fmin_cg(cost_function,
                 theta_initial,
                 fprime=cost_gradient,
                 args=args,
                 maxiter=maxiter,
                 full_output=True)
    # ============================================================
    theta_opt, cost_min, _, _, _ = ret
    print("Cost at theta found by fmin_cg: ", cost_min)
    print("theta_op: \n", theta_opt)

    # 绘制分类面
    plot_decision_boundary(theta_opt, X, y)
    plt.show()

    # 预测考试一得45分，考试二得85分的学生的录取概率
    x_test = np.array([1, 45, 85.0])
    prob = sigmoid(np.dot(theta_opt, x_test))
    print('For a student with scores 45 and 85, we predict an admission probability of: ', prob)

    # 计算在训练集上的分类正确率
    p = predict(theta_opt, X)
    print("Train Accuracy: ", np.mean(p == y)*100.)
    
logistic_regression()

######################### cell 24 #########################

# 加载数据
data = np.loadtxt("LR_data2.txt", delimiter=",")
X, y = data[:, :2], data[:, 2]
# print(X.shape) # (118, 2)
# print(y.shape) # (118,)

# 可视化数据
# ====================== YOUR CODE HERE ================
plot_data(X, y)
# ======================================================
plt.show()

######################### cell 26 #########################

def map_feature(X1, X2, degree=6):
    """Feature mapping function to polynomial features."""
    m = len(X1)
    assert len(X1) == len(X2)
    n = int((degree+2)*(degree+1)/2)

    out = np.zeros((m, n))

    idx = 0
    for i in range(degree+1):
        for j in range(i+1):
            # print i-j, j, idx
            out[:, idx] = np.power(X1, i-j)*np.power(X2, j)
            idx += 1

    return out

######################### cell 28 #########################

def cost_function_reg(theta, X, y, lmb):
    """逻辑回归的代价函数，有正则项。"""
    m = 1.0*len(y)
    J = 0
    # ====================== YOUR CODE HERE ======================
    h = sigmoid(X@theta)
    h = np.clip(h, 1e-10, 1-1e-10)
    J += np.sum(-y*np.log(h) - (1-y)*np.log(1-h)) / m
    J += lmb / (2*m) * np.sum(theta[1:]**2)
    # ============================================================

    return J


def cost_gradient_reg(theta, X, y, lmb):
    """逻辑回归的代价函数的梯度，有正则项。"""

    m = 1.0*len(y)
    grad = np.zeros_like(theta)
    # ====================== YOUR CODE HERE ======================
    h = sigmoid(X@theta)
    grad = (X.T@(h-y)) / m
    grad[1:] = grad[1:] + (lmb/m) * theta[1:]
    # ============================================================

    return grad

######################### cell 30 #########################

def logistic_regression_reg(lmb=1.0):
    """针对第二组数据建立逻辑回归模型。"""

    # 加载数据
    data = np.loadtxt("LR_data2.txt", delimiter=",")
    X, y = data[:, :2], data[:, 2]

    # 计算具有正则项的代价与梯度

    # 注意map_feature会自动加入一列 1
    X = map_feature(X[:, 0], X[:, 1])
    # print(X.shape) # (118, 28)
    
    # 初始化参数
    theta_initial = np.zeros_like(X[0, :])

    # 计算并打印初始参数对应的代价与梯度
    cost = cost_function_reg(theta_initial, X, y, lmb=lmb)
    grad = cost_gradient_reg(theta_initial, X, y, lmb=lmb)
    print("Cost at initial theta (zeros): ", cost)
    print("Gradient at initial theta (zeros): \n", grad)

    # 使用 scipy.optimize.fmin_cg 优化模型参数
    args = (X, y, lmb)
    maxiter = 200
    # ====================== YOUR CODE HERE ======================
    ret = op.fmin_cg(cost_function_reg,
                     theta_initial,
                     fprime=cost_gradient_reg,
                     args=args,
                     maxiter=maxiter,
                     full_output=True)
    # ============================================================
    theta_opt, cost_min, _, _, _ = ret
    print("Cost at theta found by fmin_cg: ", cost_min)
    print("theta_op: \n", theta_opt)

    # 绘制分类面
    plot_decision_boundary(theta_opt, X, y)
    plt.title("lambda = " + str(lmb))
    plt.show()

    # 计算在训练集上的分类正确率
    pred = predict(theta_opt, X)
    print("Train Accuracy: ", np.mean(pred == y)*100)
    
# 可选：尝试不同正则化系数lmb = 0.0, 1.0, 10.0, 100.0对分类面的影响
logistic_regression_reg(lmb=1.0) 

LR_data1.txt

shape=(100, 3)
X: (100, 2)
Y: (100,)

34.62365962451697,78.0246928153624,0
30.28671076822607,43.89499752400101,0
35.84740876993872,72.90219802708364,0
60.18259938620976,86.30855209546826,1
79.0327360507101,75.3443764369103,1
45.08327747668339,56.3163717815305,0
61.10666453684766,96.51142588489624,1
75.02474556738889,46.55401354116538,1
76.09878670226257,87.42056971926803,1
84.43281996120035,43.53339331072109,1
95.86155507093572,38.22527805795094,0
75.01365838958247,30.60326323428011,0
82.30705337399482,76.48196330235604,1
69.36458875970939,97.71869196188608,1
39.53833914367223,76.03681085115882,0
53.9710521485623,89.20735013750205,1
69.07014406283025,52.74046973016765,1
67.94685547711617,46.67857410673128,0
70.66150955499435,92.92713789364831,1
76.97878372747498,47.57596364975532,1
67.37202754570876,42.83843832029179,0
89.67677575072079,65.79936592745237,1
50.534788289883,48.85581152764205,0
34.21206097786789,44.20952859866288,0
77.9240914545704,68.9723599933059,1
62.27101367004632,69.95445795447587,1
80.1901807509566,44.82162893218353,1
93.114388797442,38.80067033713209,0
61.83020602312595,50.25610789244621,0
38.78580379679423,64.99568095539578,0
61.379289447425,72.80788731317097,1
85.40451939411645,57.05198397627122,1
52.10797973193984,63.12762376881715,0
52.04540476831827,69.43286012045222,1
40.23689373545111,71.16774802184875,0
54.63510555424817,52.21388588061123,0
33.91550010906887,98.86943574220611,0
64.17698887494485,80.90806058670817,1
74.78925295941542,41.57341522824434,0
34.1836400264419,75.2377203360134,0
83.90239366249155,56.30804621605327,1
51.54772026906181,46.85629026349976,0
94.44336776917852,65.56892160559052,1
82.36875375713919,40.61825515970618,0
51.04775177128865,45.82270145776001,0
62.22267576120188,52.06099194836679,0
77.19303492601364,70.45820000180959,1
97.77159928000232,86.7278223300282,1
62.07306379667647,96.76882412413983,1
91.56497449807442,88.69629254546599,1
79.94481794066932,74.16311935043758,1
99.2725269292572,60.99903099844988,1
90.54671411399852,43.39060180650027,1
34.52451385320009,60.39634245837173,0
50.2864961189907,49.80453881323059,0
49.58667721632031,59.80895099453265,0
97.64563396007767,68.86157272420604,1
32.57720016809309,95.59854761387875,0
74.24869136721598,69.82457122657193,1
71.79646205863379,78.45356224515052,1
75.3956114656803,85.75993667331619,1
35.28611281526193,47.02051394723416,0
56.25381749711624,39.26147251058019,0
30.05882244669796,49.59297386723685,0
44.66826172480893,66.45008614558913,0
66.56089447242954,41.09209807936973,0
40.45755098375164,97.53518548909936,1
49.07256321908844,51.88321182073966,0
80.27957401466998,92.11606081344084,1
66.74671856944039,60.99139402740988,1
32.72283304060323,43.30717306430063,0
64.0393204150601,78.03168802018232,1
72.34649422579923,96.22759296761404,1
60.45788573918959,73.09499809758037,1
58.84095621726802,75.85844831279042,1
99.82785779692128,72.36925193383885,1
47.26426910848174,88.47586499559782,1
50.45815980285988,75.80985952982456,1
60.45555629271532,42.50840943572217,0
82.22666157785568,42.71987853716458,0
88.9138964166533,69.80378889835472,1
94.83450672430196,45.69430680250754,1
67.31925746917527,66.58935317747915,1
57.23870631569862,59.51428198012956,1
80.36675600171273,90.96014789746954,1
68.46852178591112,85.59430710452014,1
42.0754545384731,78.84478600148043,0
75.47770200533905,90.42453899753964,1
78.63542434898018,96.64742716885644,1
52.34800398794107,60.76950525602592,0
94.09433112516793,77.15910509073893,1
90.44855097096364,87.50879176484702,1
55.48216114069585,35.57070347228866,0
74.49269241843041,84.84513684930135,1
89.84580670720979,45.35828361091658,1
83.48916274498238,48.38028579728175,1
42.2617008099817,87.10385094025457,1
99.31500880510394,68.77540947206617,1
55.34001756003703,64.9319380069486,1
74.77589300092767,89.52981289513276,1

LR_data2.txt

shape=(118, 3)
X: (118, 2)
Y: (118,)

0.051267,0.69956,1
-0.092742,0.68494,1
-0.21371,0.69225,1
-0.375,0.50219,1
-0.51325,0.46564,1
-0.52477,0.2098,1
-0.39804,0.034357,1
-0.30588,-0.19225,1
0.016705,-0.40424,1
0.13191,-0.51389,1
0.38537,-0.56506,1
0.52938,-0.5212,1
0.63882,-0.24342,1
0.73675,-0.18494,1
0.54666,0.48757,1
0.322,0.5826,1
0.16647,0.53874,1
-0.046659,0.81652,1
-0.17339,0.69956,1
-0.47869,0.63377,1
-0.60541,0.59722,1
-0.62846,0.33406,1
-0.59389,0.005117,1
-0.42108,-0.27266,1
-0.11578,-0.39693,1
0.20104,-0.60161,1
0.46601,-0.53582,1
0.67339,-0.53582,1
-0.13882,0.54605,1
-0.29435,0.77997,1
-0.26555,0.96272,1
-0.16187,0.8019,1
-0.17339,0.64839,1
-0.28283,0.47295,1
-0.36348,0.31213,1
-0.30012,0.027047,1
-0.23675,-0.21418,1
-0.06394,-0.18494,1
0.062788,-0.16301,1
0.22984,-0.41155,1
0.2932,-0.2288,1
0.48329,-0.18494,1
0.64459,-0.14108,1
0.46025,0.012427,1
0.6273,0.15863,1
0.57546,0.26827,1
0.72523,0.44371,1
0.22408,0.52412,1
0.44297,0.67032,1
0.322,0.69225,1
0.13767,0.57529,1
-0.0063364,0.39985,1
-0.092742,0.55336,1
-0.20795,0.35599,1
-0.20795,0.17325,1
-0.43836,0.21711,1
-0.21947,-0.016813,1
-0.13882,-0.27266,1
0.18376,0.93348,0
0.22408,0.77997,0
0.29896,0.61915,0
0.50634,0.75804,0
0.61578,0.7288,0
0.60426,0.59722,0
0.76555,0.50219,0
0.92684,0.3633,0
0.82316,0.27558,0
0.96141,0.085526,0
0.93836,0.012427,0
0.86348,-0.082602,0
0.89804,-0.20687,0
0.85196,-0.36769,0
0.82892,-0.5212,0
0.79435,-0.55775,0
0.59274,-0.7405,0
0.51786,-0.5943,0
0.46601,-0.41886,0
0.35081,-0.57968,0
0.28744,-0.76974,0
0.085829,-0.75512,0
0.14919,-0.57968,0
-0.13306,-0.4481,0
-0.40956,-0.41155,0
-0.39228,-0.25804,0
-0.74366,-0.25804,0
-0.69758,0.041667,0
-0.75518,0.2902,0
-0.69758,0.68494,0
-0.4038,0.70687,0
-0.38076,0.91886,0
-0.50749,0.90424,0
-0.54781,0.70687,0
0.10311,0.77997,0
0.057028,0.91886,0
-0.10426,0.99196,0
-0.081221,1.1089,0
0.28744,1.087,0
0.39689,0.82383,0
0.63882,0.88962,0
0.82316,0.66301,0
0.67339,0.64108,0
1.0709,0.10015,0
-0.046659,-0.57968,0
-0.23675,-0.63816,0
-0.15035,-0.36769,0
-0.49021,-0.3019,0
-0.46717,-0.13377,0
-0.28859,-0.060673,0
-0.61118,-0.067982,0
-0.66302,-0.21418,0
-0.59965,-0.41886,0
-0.72638,-0.082602,0
-0.83007,0.31213,0
-0.72062,0.53874,0
-0.59389,0.49488,0
-0.48445,0.99927,0
-0.0063364,0.99927,0
0.63265,-0.030612,0