Python数据分析与机器学习-逻辑回归案例分析

Logistic Regression

The Data

我们将建立一个逻辑回归模型来预测一个学生是否被大学录取。假设你是一个大学系的管理员，你想根据两次考试的结果来决定每个申请人的录取机会。你有以前的申请人的历史数据，你可以用它作为逻辑回归的训练集。对于每一个培训例子，你有两个考试的申请人的分数和录取决定。为了做到这一点，我们将建立一个分类模型，根据考试成绩估计入学概率。

# 三大件
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

import os
path = 'data'+os.sep+'LogiReg_data.txt'
pdData = pd.read_csv(path,header=None,names=['Exam 1','Exam 2','Admitted'])
pdData.head()

	Exam 1	Exam 2	Admitted
0	34.623660	78.024693	0
1	30.286711	43.894998	0
2	35.847409	72.902198	0
3	60.182599	86.308552	1
4	79.032736	75.344376	1

pdData.shape

(100, 3)

positive = pdData[pdData['Admitted']==1]
negative = pdData[pdData['Admitted']==0]
negative.head()

fig,ax = plt.subplots(figsize=(10,5))
ax.scatter(positive['Exam 1'],positive['Exam 2'],s=30,c='b',marker='o',label='Admitted')
ax.scatter(negative['Exam 1'],negative['Exam 2'],s=30,c='r',marker='x',label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam 1 Score')
ax.set_ylabel('Exam 2 Score')

Text(0, 0.5, 'Exam 2 Score')

The logistic regression

目标：建立分类器（求解出三个参数( heta_0, heta_1, heta_2)）

设定阈值，根据阈值判断录取结果

要完成的模块

• sigmoid：映射到概率的函数
• model：返回预测结果值
• cost：根据参数计算损失
• gradient：计算每个参数的梯度方向
• descent：进行参数更新
• accuracy：计算精度

Sigmoid 函数

[g(z) = frac{1}{1+e^{-z}} ]

• (g:mathbb{R} o [0,1])
• (g(0)=0.5)
• (g(- infty)=0)
• (g(+ infty)=1)

def sigmoid(z):
    return 1/(1+np.exp(-z))

nums = np.arange(-10,10,step=1)
print(nums)
fig,ax = plt.subplots(figsize=(12,4))
ax.plot(nums,sigmoid(nums),'r')
plt.show()

[-10  -9  -8  -7  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6   7
   8   9]

def model(X,theta):
    """ Returns our model result
    :param X: examples to classify, m x p
    :param theta: parameters,  p x 1
    :return: the sigmoid evaluated for each examples in X given parameters theta as a m x 1 vector
    """
    return sigmoid(np.matmul(X,theta))

[egin{array}{ccc} egin{pmatrix} 1 & x_{1} & x_{2}end{pmatrix} & imes & egin{pmatrix} heta_{0}\ heta_{1}\ heta_{2} end{pmatrix}end{array}= heta_{0}+ heta_{1}x_{1}+ heta_{2}x_{2} ]

pdData.insert(0,'Ones',1)
pdData.head()

# set X (training data) and y (target variable)
orig_data = pdData.as_matrix() # convert the Pandas representation of the data to an array useful for further computations
print(orig_data)
cols = orig_data.shape[1]
print(cols)
X = orig_data[:,0:cols-1]
y = orig_data[:,cols-1:cols]

# print(X[:5])
theta = np.zeros([cols-1,1])

[[ 1.         34.62365962 78.02469282  0.        ]
 [ 1.         30.28671077 43.89499752  0.        ]
 [ 1.         35.84740877 72.90219803  0.        ]
 [ 1.         60.18259939 86.3085521   1.        ]
 [ 1.         79.03273605 75.34437644  1.        ]
 [ 1.         45.08327748 56.31637178  0.        ]
 [ 1.         61.10666454 96.51142588  1.        ]
 [ 1.         75.02474557 46.55401354  1.        ]
 [ 1.         76.0987867  87.42056972  1.        ]
 [ 1.         84.43281996 43.53339331  1.        ]
 [ 1.         95.86155507 38.22527806  0.        ]
 [ 1.         75.01365839 30.60326323  0.        ]
 [ 1.         82.30705337 76.4819633   1.        ]
 [ 1.         69.36458876 97.71869196  1.        ]
 [ 1.         39.53833914 76.03681085  0.        ]
 [ 1.         53.97105215 89.20735014  1.        ]
 [ 1.         69.07014406 52.74046973  1.        ]
 [ 1.         67.94685548 46.67857411  0.        ]
 [ 1.         70.66150955 92.92713789  1.        ]
 [ 1.         76.97878373 47.57596365  1.        ]
 [ 1.         67.37202755 42.83843832  0.        ]
 [ 1.         89.67677575 65.79936593  1.        ]
 [ 1.         50.53478829 48.85581153  0.        ]
 [ 1.         34.21206098 44.2095286   0.        ]
 [ 1.         77.92409145 68.97235999  1.        ]
 [ 1.         62.27101367 69.95445795  1.        ]
 [ 1.         80.19018075 44.82162893  1.        ]
 [ 1.         93.1143888  38.80067034  0.        ]
 [ 1.         61.83020602 50.25610789  0.        ]
 [ 1.         38.7858038  64.99568096  0.        ]
 [ 1.         61.37928945 72.80788731  1.        ]
 [ 1.         85.40451939 57.05198398  1.        ]
 [ 1.         52.10797973 63.12762377  0.        ]
 [ 1.         52.04540477 69.43286012  1.        ]
 [ 1.         40.23689374 71.16774802  0.        ]
 [ 1.         54.63510555 52.21388588  0.        ]
 [ 1.         33.91550011 98.86943574  0.        ]
 [ 1.         64.17698887 80.90806059  1.        ]
 [ 1.         74.78925296 41.57341523  0.        ]
 [ 1.         34.18364003 75.23772034  0.        ]
 [ 1.         83.90239366 56.30804622  1.        ]
 [ 1.         51.54772027 46.85629026  0.        ]
 [ 1.         94.44336777 65.56892161  1.        ]
 [ 1.         82.36875376 40.61825516  0.        ]
 [ 1.         51.04775177 45.82270146  0.        ]
 [ 1.         62.22267576 52.06099195  0.        ]
 [ 1.         77.19303493 70.4582      1.        ]
 [ 1.         97.77159928 86.72782233  1.        ]
 [ 1.         62.0730638  96.76882412  1.        ]
 [ 1.         91.5649745  88.69629255  1.        ]
 [ 1.         79.94481794 74.16311935  1.        ]
 [ 1.         99.27252693 60.999031    1.        ]
 [ 1.         90.54671411 43.39060181  1.        ]
 [ 1.         34.52451385 60.39634246  0.        ]
 [ 1.         50.28649612 49.80453881  0.        ]
 [ 1.         49.58667722 59.80895099  0.        ]
 [ 1.         97.64563396 68.86157272  1.        ]
 [ 1.         32.57720017 95.59854761  0.        ]
 [ 1.         74.24869137 69.82457123  1.        ]
 [ 1.         71.79646206 78.45356225  1.        ]
 [ 1.         75.39561147 85.75993667  1.        ]
 [ 1.         35.28611282 47.02051395  0.        ]
 [ 1.         56.2538175  39.26147251  0.        ]
 [ 1.         30.05882245 49.59297387  0.        ]
 [ 1.         44.66826172 66.45008615  0.        ]
 [ 1.         66.56089447 41.09209808  0.        ]
 [ 1.         40.45755098 97.53518549  1.        ]
 [ 1.         49.07256322 51.88321182  0.        ]
 [ 1.         80.27957401 92.11606081  1.        ]
 [ 1.         66.74671857 60.99139403  1.        ]
 [ 1.         32.72283304 43.30717306  0.        ]
 [ 1.         64.03932042 78.03168802  1.        ]
 [ 1.         72.34649423 96.22759297  1.        ]
 [ 1.         60.45788574 73.0949981   1.        ]
 [ 1.         58.84095622 75.85844831  1.        ]
 [ 1.         99.8278578  72.36925193  1.        ]
 [ 1.         47.26426911 88.475865    1.        ]
 [ 1.         50.4581598  75.80985953  1.        ]
 [ 1.         60.45555629 42.50840944  0.        ]
 [ 1.         82.22666158 42.71987854  0.        ]
 [ 1.         88.91389642 69.8037889   1.        ]
 [ 1.         94.83450672 45.6943068   1.        ]
 [ 1.         67.31925747 66.58935318  1.        ]
 [ 1.         57.23870632 59.51428198  1.        ]
 [ 1.         80.366756   90.9601479   1.        ]
 [ 1.         68.46852179 85.5943071   1.        ]
 [ 1.         42.07545454 78.844786    0.        ]
 [ 1.         75.47770201 90.424539    1.        ]
 [ 1.         78.63542435 96.64742717  1.        ]
 [ 1.         52.34800399 60.76950526  0.        ]
 [ 1.         94.09433113 77.15910509  1.        ]
 [ 1.         90.44855097 87.50879176  1.        ]
 [ 1.         55.48216114 35.57070347  0.        ]
 [ 1.         74.49269242 84.84513685  1.        ]
 [ 1.         89.84580671 45.35828361  1.        ]
 [ 1.         83.48916274 48.3802858   1.        ]
 [ 1.         42.26170081 87.10385094  1.        ]
 [ 1.         99.31500881 68.77540947  1.        ]
 [ 1.         55.34001756 64.93193801  1.        ]
 [ 1.         74.775893   89.5298129   1.        ]]
4


C:MyProgramsAnaconda3libsite-packagesipykernel_launcher.py:5: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.
  """

X[:5]

array([[ 1.        , 34.62365962, 78.02469282],
       [ 1.        , 30.28671077, 43.89499752],
       [ 1.        , 35.84740877, 72.90219803],
       [ 1.        , 60.18259939, 86.3085521 ],
       [ 1.        , 79.03273605, 75.34437644]])

y[:5]

array([[0.],
       [0.],
       [0.],
       [1.],
       [1.]])

theta

array([[0.],
       [0.],
       [0.]])

X.shape,y.shape,theta.shape

((100, 3), (100, 1), (3, 1))

损失函数

将对数似然函数去负号

[D(h_ heta(x), y) = -ylog(h_ heta(x)) - (1-y)log(1-h_ heta(x)) ]

求平均损失

[J( heta)=frac{1}{m}sum_{i=1}^{m} D(h_ heta(x_i), y_i) ]

def costFunction(X,y,theta):
    left = np.multiply(-y,np.log(model(X,theta))) # 同*，元素级乘法
    right = np.multiply((1-y),np.log(1-model(X,theta)))
    return np.sum(left-right)/(len(X))

costFunction(X,y,theta)

0.6931471805599453

计算梯度

[frac{partial J}{partial heta_j}=-frac{1}{m}sum_{i=1}^m (y_i - h_ heta (x_i))x^i_{j} ]

[egin{pmatrix}frac{partial J}{partial heta_0}\ frac{partial J}{partial heta_1}\ frac{partial J}{partial heta_2} end{pmatrix}=-frac{1}{m}egin{pmatrix} 1 & cdots & 1\ x^1_{1} & cdots & x^m_{1}\ x^1_{2}& cdots & x^m_{2}end{pmatrix}egin{pmatrix}y_1 - h_ heta (x_1)\ vdots\ y_m - h_ heta (x_m) end{pmatrix}=frac{1}{m}X^T(h_ heta(x)-y) ]

def gradient(X,y,theta):
    grad = np.zeros(theta.shape)
    # error = (model(X, theta)- y).ravel()
    # for j in range(len(theta.ravel())): #for each parmeter
    #     term = np.multiply(error, X[:,j])
    #     grad[0, j] = np.sum(term) / len(X)
    error = np.matmul(X.T,(model(X,theta)-y))
    grad = error/len(X)
    return grad

gradient(X,y,theta)

array([[ -0.1       ],
       [-12.00921659],
       [-11.26284221]])

Gradient descent

比较3种不同梯度下降方法

STOP_ITER = 0
STOP_COST = 1
STOP_GRAD = 2

def stopCriterion(type,value,threshold):
    # 设定三种不同的停止策略
    if type == STOP_ITER: return value > threshold
    elif type == STOP_COST: return abs(value[-1]-value[-2]) < threshold
    elif type == STOP_GRAD: return np.linalg.norm(value) < threshold

import time

def descent(data,theta,batchSize, stopType, thresh,alpha):
    #梯度下降求解
    
    init_time = time.time()
    i = 0 #  迭代次数
    k = 0 # batch
    X,y = shuffleData(data)
    grad = np.zeros(theta.shape) # 计算梯度
    costs = [costFunction(X,y,theta)] # 损失值
    
    while True:
        grad = gradient(X[k:k+batchSize],y[k:k+batchSize],theta)
        k += batchSize # 取batch数量个数据
        if k >= n:
            k = 0;
            X,y = shuffleData(data) # 重新洗牌
        theta = theta - alpha*grad # 参数更新
        costs.append(costFunction(X,y,theta)) # 计算新的损失
        i += 1
        
        if stopType == STOP_ITER: value = i
        elif stopType == STOP_COST: value = costs
        elif stopType == STOP_GRAD: value = grad
        if stopCriterion(stopType,value,thresh): break
        
    return theta,i-1,costs,grad,time.time()-init_time

import numpy.random
# 洗牌
def shuffleData(data):
    np.random.shuffle(data)
    cols = data.shape[1]
    X = data[:, 0:cols-1]
    y = data[:, cols-1:]
    return X, y

def runExpe(data, theta, batchSize, stopType, thresh, alpha):
    #import pdb; pdb.set_trace();
    theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh, alpha)
    name = "Original" if (data[:,1]>2).sum() > 1 else "Scaled"
    name += " data - learning rate: {} - ".format(alpha)
    if batchSize==n: strDescType = "Gradient"
    elif batchSize==1:  strDescType = "Stochastic"
    else: strDescType = "Mini-batch ({})".format(batchSize)
    name += strDescType + " descent - Stop: "
    if stopType == STOP_ITER: strStop = "{} iterations".format(thresh)
    elif stopType == STOP_COST: strStop = "costs change < {}".format(thresh)
    else: strStop = "gradient norm < {}".format(thresh)
    name += strStop
    print ("***{}
Theta: {} - Iter: {} - Last cost: {:03.2f} - Duration: {:03.2f}s".format(
        name, theta, iter, costs[-1], dur))
    fig, ax = plt.subplots(figsize=(12,4))
    ax.plot(np.arange(len(costs)), costs, 'r')
    ax.set_xlabel('Iterations')
    ax.set_ylabel('Cost')
    ax.set_title(name.upper() + ' - Error vs. Iteration')
    return theta

不同的停止策略

设定迭代次数

#选择的梯度下降方法是基于所有样本的
n=100
runExpe(orig_data, theta, n, STOP_ITER, thresh=5000, alpha=0.000001)

***Original data - learning rate: 1e-06 - Gradient descent - Stop: 5000 iterations
Theta: [[-0.00027127]
 [ 0.00705232]
 [ 0.00376711]] - Iter: 5000 - Last cost: 0.63 - Duration: 1.65s





array([[-0.00027127],
       [ 0.00705232],
       [ 0.00376711]])

根据损失值停止

设定阈值为1E-6

runExpe(orig_data, theta, n, STOP_COST, thresh=0.000001, alpha=0.001)

***Original data - learning rate: 0.001 - Gradient descent - Stop: costs change < 1e-06
Theta: [[-5.13364014]
 [ 0.04771429]
 [ 0.04072397]] - Iter: 109901 - Last cost: 0.38 - Duration: 48.19s





array([[-5.13364014],
       [ 0.04771429],
       [ 0.04072397]])

对比不同的梯度下降法

Stochastic descent

runExpe(orig_data, theta, 1, STOP_ITER, thresh=5000, alpha=0.001)

***Original data - learning rate: 0.001 - Stochastic descent - Stop: 5000 iterations
Theta: [[-0.39554857]
 [ 0.02237485]
 [-0.0276303 ]] - Iter: 5000 - Last cost: 0.89 - Duration: 0.76s





array([[-0.39554857],
       [ 0.02237485],
       [-0.0276303 ]])

有点爆炸。。。很不稳定,再来试试把学习率调小一些

runExpe(orig_data, theta, 1, STOP_ITER, thresh=15000, alpha=0.000002)

***Original data - learning rate: 2e-06 - Stochastic descent - Stop: 15000 iterations
Theta: [[-0.00202176]
 [ 0.01001261]
 [ 0.00091419]] - Iter: 15000 - Last cost: 0.63 - Duration: 1.86s





array([[-0.00202176],
       [ 0.01001261],
       [ 0.00091419]])

速度快，但稳定性差，需要很小的学习率

Mini-batch descent

runExpe(orig_data, theta, 16, STOP_ITER, thresh=15000, alpha=0.001)

***Original data - learning rate: 0.001 - Mini-batch (16) descent - Stop: 15000 iterations
Theta: [[-1.03466277]
 [ 0.02768388]
 [ 0.01858675]] - Iter: 15000 - Last cost: 0.72 - Duration: 2.89s





array([[-1.03466277],
       [ 0.02768388],
       [ 0.01858675]])

浮动仍然比较大，我们来尝试下对数据进行标准化 将数据按其属性(按列进行)减去其均值，然后除以其方差。最后得到的结果是，对每个属性/每列来说所有数据都聚集在0附近，方差值为1

from sklearn import preprocessing as pp

scaled_data = orig_data.copy()
scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3])

runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.001)

***Scaled data - learning rate: 0.001 - Gradient descent - Stop: 5000 iterations
Theta: [[0.3080807 ]
 [0.86494967]
 [0.77367651]] - Iter: 5000 - Last cost: 0.38 - Duration: 2.30s





array([[0.3080807 ],
       [0.86494967],
       [0.77367651]])

它好多了！原始数据，只能达到达到0.61，而我们得到了0.38个在这里！ 所以对数据做预处理是非常重要的

runExpe(scaled_data, theta, n, STOP_GRAD, thresh=0.02, alpha=0.001)

***Scaled data - learning rate: 0.001 - Gradient descent - Stop: gradient norm < 0.02
Theta: [[1.0707921 ]
 [2.63030842]
 [2.41079787]] - Iter: 59422 - Last cost: 0.22 - Duration: 31.39s





array([[1.0707921 ],
       [2.63030842],
       [2.41079787]])

更多的迭代次数会使得损失下降的更多！

theta = runExpe(scaled_data, theta, 1, STOP_GRAD, thresh=0.002/5, alpha=0.001)

***Scaled data - learning rate: 0.001 - Stochastic descent - Stop: gradient norm < 0.0004
Theta: [[1.14837004]
 [2.7936333 ]
 [2.5646749 ]] - Iter: 72578 - Last cost: 0.22 - Duration: 13.34s

随机梯度下降更快，但是我们需要迭代的次数也需要更多，所以还是用batch的比较合适！！！

runExpe(scaled_data, theta, 16, STOP_GRAD, thresh=0.002*2, alpha=0.001)

***Scaled data - learning rate: 0.001 - Mini-batch (16) descent - Stop: gradient norm < 0.004
Theta: [[1.17941228]
 [2.8498756 ]
 [2.62631998]] - Iter: 5690 - Last cost: 0.21 - Duration: 1.33s





array([[1.17941228],
       [2.8498756 ],
       [2.62631998]])

精度

#设定阈值
def predict(X, theta):
    return [1 if x >= 0.5 else 0 for x in model(X, theta)]

scaled_X = scaled_data[:, :3]
y = scaled_data[:, 3]
predictions = predict(scaled_X, theta)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
accuracy = (sum(map(int, correct)) % len(correct))
print ('accuracy = {0}%'.format(accuracy))

accuracy = 89%