吴恩达机器学习作业2（逻辑回归）

一.知识回顾

后续有代价函数误差什么什么的

二.题目

题目概述：在训练的初始阶段，我们将要构建一个逻辑回归模型来预测，某个学生是否被大学录取。 设想你是大学相关部分的管理者，想通过申请学生两次测试的评分，来决定他们是否被录取。 现在你拥有之前申请学生的可以用于训练逻辑回归的训练样本集。对于每一个训练样本，你有他们两次测试的评分和最后是被录取的结果。

二.编程的设计

【1】引入一些包

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('fivethirtyeight')
import matplotlib.pyplot as plt
# import tensorflow as tf
from sklearn.metrics import classification_report#这个包是评价报告

【2】准备数据

data = pd.read_csv('ex2data1.txt', names=['exam1', 'exam2', 'admitted'])
data.head()#看前五行

0无法录取，1可以录取

data.describe()

和上一篇一样，没啥好讲的

sns.set(context="notebook", style="darkgrid", palette=sns.color_palette("RdBu", 2))

sns.lmplot('exam1', 'exam2', hue='admitted', data=data, 
           size=6, 
           fit_reg=False, 
           scatter_kws={"s": 50}
          )
plt.show()#看下数据的样子

hue函数是指将name中的数据都集中显示在一张图上

def get_X(df):#读取特征
#     """
#     use concat to add intersect feature to avoid side effect
#     not efficient for big dataset though
#     """
    ones = pd.DataFrame({'ones': np.ones(len(df))})#ones是m行1列的dataframe
    data = pd.concat([ones, df], axis=1)  # 合并数据，根据列合并，concat拼接函数 ，相当于特征进行扩展
    return data.iloc[:, :-1].as_matrix()  # 这个操作返回 ndarray,不是矩阵，-1表示最后一列，但是左开右闭，所以正好是扩展后的前三列


def get_y(df):#读取标签（去最后一列）
#     '''assume the last column is the target'''
    return np.array(df.iloc[:, -1]) #df.iloc[:, -1]是指df的最后一列


def normalize_feature(df):
#     """Applies function along input axis(default 0) of DataFrame."""
    return df.apply(lambda column: (column - column.mean()) / column.std())#特征缩放 （每一列的元素-每一列的均值）/每一列的方差

X = get_X(data)
print(X.shape)

y = get_y(data)
print(y.shape)

(100, 3)
(100,) # 这是一个列向量

【3】sigmoid 函数

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为： $g\left(z\right)=\frac{1}{1+e^{-z}}$
合起来，我们得到逻辑回归模型的假设函数：
$h_{\theta }\left(x\right)=\frac{1}{1+e^{-{\theta }^{T}X}}$

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

fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(np.arange(-10, 10, step=0.01),
        sigmoid(np.arange(-10, 10, step=0.01)))
ax.set_ylim((-0.1,1.1))
ax.set_xlabel('z', fontsize=18)
ax.set_ylabel('g(z)', fontsize=18)
ax.set_title('sigmoid function', fontsize=18)
plt.show()

【4】cost function(代价函数)

• $max(\ell (\theta ))=min(-\ell (\theta ))$
• choose $-\ell (\theta )$ as the cost function

1.设置参数

theta = theta=np.zeros(3) # X(m*n) so theta is n*1
theta

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

def cost(theta, X, y):
    ''' cost fn is -l(theta) for you to minimize'''
    return np.mean(-y * np.log(sigmoid(X @ theta)) - (1 - y) * np.log(1 - sigmoid(X @ theta)))

# X @ theta与X.dot(theta)等价

上面这个代码的流程图：

测验打分：

cost(theta, X, y) 

# 0.69314718055994529

【5】gradient descent(梯度下降)

1.损失梯度的推导

• 这是批量梯度下降（batch gradient descent）
• 转化为向量化计算： $\frac{1}{m}{X}^T(Sigmoid(X\theta )-y)$
  $$\frac{\partial J\left(\theta \right)}{\partial {\theta }_j}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }\left(x^{(i)}\right)-y^{(i)})x_{{}_j}^{(i)}$$
  

2.代码

def gradient(theta, X, y):
#     '''just 1 batch gradient'''
    return (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y)

gradient(theta, X, y)

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

【6】拟合参数

这里我使用 scipy.optimize.minimize 去寻找参数

import scipy.optimize as opt

res = opt.minimize(fun=cost, x0=theta, args=(X, y), method='Newton-CG', jac=gradient)

print(res)

x=（参数1，参数2，参数3）这是优化后的函数

【7】用训练集预测和验证

def predict(x, theta):
    prob = sigmoid(x @ theta)
    return (prob >= 0.5).astype(int)

final_theta = res.x
y_pred = predict(X, final_theta)

print(classification_report(y, y_pred))

查准率（Precision）和查全率（Recall）
见：
这两个越高越好
具体介绍见我的笔记吴恩达机器学习day8（机器学习系统的设计）

【7】寻找决策边界

$X\times \theta =0$ (this is the line)

1.输出最后拟合的参数

print(res.x) # this is final theta

# [-25.16227358 0.20623923 0.20147921]

2.根据边界公式计算$X\times \theta$

coef = -(res.x / res.x[2])  # find the equation
print(coef)

x = np.arange(130, step=0.1)
y = coef[0] + coef[1]*x

# [ 124.88769463 -1.0236254 -1. ]

data.describe()  # find the range of x and y

3.画出图形

sns.set(context="notebook", style="ticks", font_scale=1.5)

sns.lmplot('exam1', 'exam2', hue='admitted', data=data, 
           size=6, 
           fit_reg=False, 
           scatter_kws={"s": 25}
          )

plt.plot(x, y, 'grey')
plt.xlim(0, 130)
plt.ylim(0, 130)
plt.title('Decision Boundary')
plt.show()

刚好可以用一条直线分开

三. 正则化逻辑回归

df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])
df.head()

sns.set(context="notebook", style="ticks", font_scale=1.5)

sns.lmplot('test1', 'test2', hue='accepted', data=df, 
           size=6, 
           fit_reg=False, 
           scatter_kws={"s": 50}
          )

plt.title('Regularized Logistic Regression')
plt.show()

由于无法用一个直线区分，我们是无法用逻辑回归分类哦，所以要用特征映射，将低微映射映射到高维特征

特征映射

对特征进行扩张

x1 = np.array(df.test1)
x2 = np.array(df.test2)

data = feature_mapping(x1, x2, power=6)
print(data.shape)
data.head() #(118, 28)

正则化函数

$$J\left(\theta \right)=\frac{1}{m}\sum _{i=1}^m[-y^{(i)}\log \left(h_{\theta }\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right)\log \left(1-h_{\theta }\left(x^{(i)}\right)\right)]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$

theta = np.zeros(data.shape[1])
X = feature_mapping(x1, x2, power=6, as_ndarray=True)
print(X.shape)

y = get_y(df)
print(y.shape)

(118, 28)
(118,)

def regularized_cost(theta, X, y, l=1):
#     '''you don't penalize theta_0'''
    theta_j1_to_n = theta[1:]
    regularized_term = (l / (2 * len(X))) * np.power(theta_j1_to_n, 2).sum()

    return cost(theta, X, y) + regularized_term
#正则化代价函数

regularized_cost(theta, X, y, l=1)

0.6931471805599454
因为我们设置theta为0，所以这个正则化代价函数与代价函数的值相同

regularized gradient(正则化梯度)

$$\frac{\partial J\left(\theta \right)}{\partial {\theta }_j}=\left(\frac{1}{m}\sum _{i=1}^m\left(h_{\theta }\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)\right)+\frac{\lambda }{m}{\theta }_j\text{ for j}\ge \text{1}$$

def regularized_gradient(theta, X, y, l=1):
#     '''still, leave theta_0 alone'''
    theta_j1_to_n = theta[1:]
    regularized_theta = (l / len(X)) * theta_j1_to_n

    # by doing this, no offset is on theta_0
    regularized_term = np.concatenate([np.array([0]), regularized_theta])

    return gradient(theta, X, y) + regularized_term

regularized_gradient(theta, X, y)

拟合参数

import scipy.optimize as opt

print('init cost = {}'.format(regularized_cost(theta, X, y)))

res = opt.minimize(fun=regularized_cost, x0=theta, args=(X, y), method='Newton-CG', jac=regularized_gradient)
res

预测

final_theta = res.x
y_pred = predict(X, final_theta)

print(classification_report(y, y_pred))

使用不同的 𝜆 （这个是常数）

画出决策边界

我们找到所有满足 𝑋×𝜃=0 的x
instead of solving polynomial equation, just create a coridate x,y grid that is dense enough, and find all those 𝑋×𝜃 that is close enough to 0, then plot them

def draw_boundary(power, l):
#     """
#     power: polynomial power for mapped feature
#     l: lambda constant
#     """
    density = 1000
    threshhold = 2 * 10**-3

    final_theta = feature_mapped_logistic_regression(power, l)
    x, y = find_decision_boundary(density, power, final_theta, threshhold)

    df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])
    sns.lmplot('test1', 'test2', hue='accepted', data=df, size=6, fit_reg=False, scatter_kws={"s": 100})

    plt.scatter(x, y, c='R', s=10)
    plt.title('Decision boundary')
    plt.show()

def feature_mapped_logistic_regression(power, l):
#     """for drawing purpose only.. not a well generealize logistic regression
#     power: int
#         raise x1, x2 to polynomial power
#     l: int
#         lambda constant for regularization term
#     """
    df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])
    x1 = np.array(df.test1)
    x2 = np.array(df.test2)
    y = get_y(df)

    X = feature_mapping(x1, x2, power, as_ndarray=True)
    theta = np.zeros(X.shape[1])

    res = opt.minimize(fun=regularized_cost,
                       x0=theta,
                       args=(X, y, l),
                       method='TNC',
                       jac=regularized_gradient)
    final_theta = res.x

    return final_theta

def find_decision_boundary(density, power, theta, threshhold):
    t1 = np.linspace(-1, 1.5, density)
    t2 = np.linspace(-1, 1.5, density)

    cordinates = [(x, y) for x in t1 for y in t2]
    x_cord, y_cord = zip(*cordinates)
    mapped_cord = feature_mapping(x_cord, y_cord, power)  # this is a dataframe

    inner_product = mapped_cord.as_matrix() @ theta

    decision = mapped_cord[np.abs(inner_product) < threshhold]

    return decision.f10, decision.f01
#寻找决策边界函数

draw_boundary(power=6, l=1)#lambda=1

draw_boundary(power=6, l=0)  # no regularization, over fitting，#lambda=0,没有正则化，过拟合了

draw_boundary(power=6, l=100)  # underfitting，#lambda=100,欠拟合