机器学习基础算法在线练习——线性回归&逻辑回归

1、线性回归

1.1 了解线性回归

线性回归是利用数理统计中回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法，运用十分广泛。其表达形式为y = w’x+e，e为误差服从均值为0的正态分布。
回归分析中，只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示，这种回归分析称为一元线性回归分析。如果回归分析中包括两个或两个以上的自变量，且因变量和自变量之间是线性关系，则称为多元线性回归分析。

1.2 背景

假设你是一家餐厅的领导，正在考虑在不同的城市开设新的分店。该连锁店已经在不同的城市有了餐车，并且你能够获得每个城市的人口和利润数据。

现在需要使用这些数据来帮助你选择下一个被扩展的城市。

文件ex1data1.txt包含线性回归问题的数据集。第一列数据对应城市人口，第二列数据对应那座城市的餐车的利润。利润为负时表示亏损。

1.3 代码&部分思路

在开始之前，对数据进行可视化通常很有用。对于该数据集，可以使用散点图进行可视化，因为它只有两个属性（人口、利润）。

1.3.1 数据可视化

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os

# 数据存储路径
path = 'ex1data1.txt'

# 读入相应的数据文件
data = pd.read_csv(filepath_or_buffer=path, header=None,names=['Population','Profit'])

###在这里填入代码###
data.plot(kind='scatter', x='Population', y='Profit',c='blue', figsize=(12,8))

图片：

1.3.2 梯度下降

我们将使用梯度下降来选择合适的线性回归参数θ用以拟合给定数据集。

线性回归的目的是最小化成本函数：

假设$h_{\theta }(X)$由以下线性模型给出：

模型的参数是${\theta }_j$的值，这些将用来调整以最小化成本$J(\theta )$。

其中一种方法是使用批量梯度下降算法，在批量梯度下降中，每次迭代地执行更新，随着梯度下降的每一步计算，参数${\theta }_j$越来越接近能够使得成本$J(\theta )$达到最低的最佳值。

(同时更新所有的${\theta }_j$）

接下来，我们在数据中添加了一个维度来拟合截距项${\theta }_0$。并将初始参数值设为0，学习率$\alpha$设为0.01。

#在列索引为0处添加数据列，该列值均为1
data.insert(0, 'Ones', 1)

#获取数据列数
cols = data.shape[1]

#对变量X和y进行初始化,并将其数据类型转换为矩阵
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1:cols]
X = np.matrix(X.values)
y = np.matrix(y.values)

#学习率、迭代次数的初始化
alpha = 0.01
iterations = 1500

1.3.3 计算成本J(θ)

在执行梯度下降最小化成本函数$J(\theta )$时，通过计算成本来监视收敛状态是有帮助的。

实现一个计算成本$J(\theta )$的函数computeCost，用于检查梯度下降实现的收敛性。

其中，X和y不是标量值，而是矩阵，其行代表训练集中的示例。

def computeCost(X, y, theta):
    inner = np.power(((X * theta.T) - y), 2)
    return np.sum(inner) / (2 * len(X))

theta = np.matrix(np.array([0,0]))
computeCost(X, y, theta)

1.3.4 梯度下降

接下来，将实现梯度下降，给出的代码已经实现了循环结构，只需要在每次的迭代中提供$\theta$的更新。

请记住，成本$J(\theta )$为参数-被向量$\theta$终止，而不是$X$和$y$。也就是说，我们将$J(\theta )$的值最小化通过改变矢量$\theta$的值，而不是通过改变$X$或$y$。

验证梯度下降是否正常工作的一种好方法是查看$J(\theta )$的值，并检查该值每步是否减小。每次迭代时，代码都会调用computeCost函数并打印成本。假设你实现了梯度下降，正确地计算成本，$J(\theta )$值永远不会增加，并且应该在算法结束时收敛到稳定值。

def gradientDescent(X, y, theta, alpha, iters):
    temp = np.matrix(np.zeros(theta.shape))
    parameters = int(theta.ravel().shape[1])
    cost = np.zeros(iters)
    
    for i in range(iters):
        error = (X * theta.T) - y
        
        for j in range(parameters):
            term = np.multiply(error, X[:,j])
            temp[0,j] = theta[0,j] - ((alpha / len(X)) * np.sum(term))
            
        theta = temp
        cost[i] = computeCost(X, y, theta)
        
    return theta, cost

g, cost = gradientDescent(X, y, theta, alpha, iterations)

#对拟合曲线进行绘制
x = np.linspace(data.Population.min(), data.Population.max(), 100)
f = g[0, 0] + (g[0, 1] * x)

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x, f, 'b', label='Prediction')
ax.scatter(data.Population, data.Profit, c='red',label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')

1.3.5 可视化成本函数

为了更好地理解成本函数的迭代计算，将每一步计算的cost值进行记录并绘制。

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iterations), cost, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')

Text(0.5, 1.0, ‘Error vs. Training Epoch’)

2、逻辑回归

2.1 了解逻辑回归：

logistic回归又称logistic回归分析，是一种广义的线性回归分析模型，常用于数据挖掘，疾病自动诊断，经济预测等领域。例如，探讨引发疾病的危险因素，并根据危险因素预测疾病发生的概率等。以胃癌病情分析为例，选择两组人群，一组是胃癌组，一组是非胃癌组，两组人群必定具有不同的体征与生活方式等。因此因变量就为是否胃癌，值为“是”或“否”，自变量就可以包括很多了，如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的，也可以是分类的。然后通过logistic回归分析，可以得到自变量的权重，从而可以大致了解到底哪些因素是胃癌的危险因素。同时根据该权值可以根据危险因素预测一个人患癌症的可能性。

2.2 背景：

假设你是大学某个部门的负责人，你要根据两次考试的结果来决定每个申请人的入学机会。目前已经有了以往申请者的历史数据，并且可以用作逻辑回归的训练集。对于每行数据，都包含对应申请者的两次考试分数和最终的录取结果。

在这里，你需要建立一个分类模型，根据这两次的考试分数来预测申请者的录取结果。

2.3 代码&部分思路：

2.3.1 数据可视化

和上面一样，先进行数据的可视化

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path = 'ex2data1.txt'
data = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])

positive = data[data['Admitted'].isin([1])]
negative = data[data['Admitted'].isin([0])]

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

2.3.2 Sigmoid函数

在正式开始之前，我们先来了解一个函数：Sigmoid函数。
我们还记得逻辑回归假设的定义是：

其中 g 代表一个常用的逻辑函数为S形函数（Sigmoid function），公式为：

合起来，我们得到逻辑回归模型的假设函数：

接下来，我们就需要编写代码实现Sigmoid函数，编写后试着测试一些值，如果x的正值较大，则函数值应接近1；如果x的负值较大，则函数值应接近0。而对于x等于0时，则函数值为0.5。

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

nums = np.arange(-10, 10, step=1)

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(nums, sigmoid(nums), 'r')
plt.show()

此时运行代码的结果应该为：

2.3.3 代价函数和梯度

①、代价函数

我们知道逻辑回归的代价函数是：

现在，我们编写代码实现代价函数以进行逻辑回归的成本计算，并且经过所给数据测试后，初始的成本约为0.693。

###在这里填入代码###
def cost(theta,X,y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y,np.log(sigmoid(X* theta.T)))
    second = np.multiply((1-y),np.log(1-sigmoid(X*theta.T)))
    return np.sum(first-second)/(len(X))

#增加一列值为1，这和我们在练习1中的操作很相似
data.insert(0, 'Ones', 1)

# 定义X为训练数据，y为目的变量
cols = data.shape[1]  # data为前面从文件中获取到的数据
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1:cols]

# 将X,y转换为numpy数组，并初始化theta值为0
X = np.array(X.values)
y = np.array(y.values)
theta = np.zeros(3)

cost(theta, X, y)

②、梯度下降

接下来，我们需要编写代码实现梯度下降用来计算我们的训练数据、标签和一些参数$\theta$的梯度。

批量梯度下降转化为向量化计算： $\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)}$$
这里需要注意的是，我们实际上没有在这个函数中执行梯度下降，我们仅仅在计算一个梯度步长。由于我们使用Python，我们可以用SciPy的optimize命名空间来做同样的事情。

def gradient(theta, X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        grad[i] = np.sum(term) / len(X)
    
    return grad

2.3.4 寻找最优参数

现在可以用SciPy’s truncated newton（TNC）实现寻找最优参数。

import scipy.optimize as opt
result = opt.fmin_tnc(func=cost, x0=theta, fprime=gradient, args=(X, y))

在这里，我们可以得到result的结果：(array([-25.16131868, 0.20623159, 0.20147149]), 36, 0)

2.3.5 评估逻辑回归

接下来，我们需要编写代码实现预测函数，用所学的最优参数$\theta$来为数据集X输出预测结果。然后，可以使用这个函数来给我们定义的分类器的训练精度进行打分。

逻辑回归的假设函数：

当$h_{\theta }$大于等于0.5时，预测 y=1

当$h_{\theta }$小于0.5时，预测 y=0。

###在这里填入代码###
def predict(theta, X):
    probability = sigmoid(X * theta.T)
    return [1 if x >= 0.5 else 0 for x in probability]

predictions = predict(theta_min, X)
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))

最后的结果应该是89%

2.4 正则化逻辑回归

在本部分练习中，我们将要通过加入正则项提升逻辑回归算法。

正则化是成本函数中的一个术语，它使算法更倾向于“更简单”的模型。这个理论助于减少过拟合，提高模型的泛化能力。

设想你是工厂的生产主管，你有一些芯片在两次测试中的测试结果。对于这两次测试，你想决定芯片是要被接受或抛弃。为了帮助你做出艰难的决定，你拥有过去芯片的测试数据集，从其中你可以构建一个逻辑回归模型。

2.4.1 数据可视化

与第一部分的练习类似，首先对数据进行可视化：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path =  'ex2data2.txt'
data2 = pd.read_csv(path, header=None, names=['Test 1', 'Test 2', 'Accepted'])

positive = data2[data2['Accepted'].isin([1])]
negative = data2[data2['Accepted'].isin([0])]

fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['Test 1'], positive['Test 2'], s=50, c='b', marker='o', label='Accepted')
ax.scatter(negative['Test 1'], negative['Test 2'], s=50, c='r', marker='x', label='Rejected')
ax.legend()
ax.set_xlabel('Test 1 Score')
ax.set_ylabel('Test 2 Score')
plt.show()

对于这部分数据，我们可以看出不同类别的数据点之间没有明显的线性决策界限用于划分两类数据。

因此，逻辑回归无法在此数据集上得到较好的效果，因为逻辑回归只能知道线性决策边界。

2.4.2 特征映射

一种能够更好地拟合数据的方法是构造从原始特征的多项式中得到的特征，即特征映射。如下图所示，作为这种映射的结果，我们的两个特征向量$x_1,x_2$（两次质量保证测试的分数）已经被转换成了28维的向量。

在这个高维特征向量上训练的逻辑回归分类器将具有更复杂的决策边界，并在二维图中绘制时呈现非线性的划分曲线。

虽然特征映射允许我们构建一个更具有表现力的分类器，但它也更容易过拟合。接下来，需要实现正则化逻辑回归用于拟合数据，并使用正则化来帮助解决过拟合问题。

我们通过创建一组多项式特征来开始！

# 设定映射深度
degree = 5
# 分别取两次测试的分数
x1 = data2['Test 1']
x2 = data2['Test 2']

data2.insert(3, 'Ones', 1)

# 设定计算方式进行映射
for i in range(1, degree):
    for j in range(0, i):
        data2['F' + str(i) + str(j)] = np.power(x1, i-j) * np.power(x2, j)

# 整理数据列
data2.drop('Test 1', axis=1, inplace=True)
data2.drop('Test 2', axis=1, inplace=True)

print("特征映射后具有特征维数：%d" %data2.shape[1])
data2.head()

特征映射后具有特征维数：12

2.4.3 代价函数和梯度

接下来，你需要编写代码来实现计算正则化逻辑回归的代价函数和梯度，并返回计算的代价值和梯度。

正则化逻辑回归的代价函数如下：
$$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$$

其中$\lambda$是“学习率”参数，其值会影响函数中的正则项值。且不应该正则化参数${\theta }_0$。

###在这里填入代码###
def costReg(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
    return np.sum(first - second) / len(X) + reg

接下来，我们需要实现正则化梯度函数，使用梯度下降法使得代价函数最小化。

因为在代价函数的计算中我们未对${\theta }_0$进行正则化，所以梯度下降算法将分为两种情况：

对上面的算法中 j=1,2,…,n 时的更新式子进行调整可得：
${\theta }_j:={\theta }_j(1-a\frac{\lambda }{m})-a\frac{1}{m}\sum _{i=1}^m(h_{\theta }\left(x^{(i)}\right)-y^{(i)})x_j^{(i)}$

###在这里填入代码###
def gradientReg(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        
        if (i == 0):
            grad[i] = np.sum(term) / len(X)
        else:
            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])
    
    return grad
    

接下来，类似于第一部分的练习中，进行变量的初始化。

# 从数据集中取得对应的特征列和标签列
cols = data2.shape[1]
X2 = data2.iloc[:,1:cols]
y2 = data2.iloc[:,0:1]

# 转换为Numpy数组并初始化theta为零矩阵
X2 = np.array(X2.values)
y2 = np.array(y2.values)
theta2 = np.zeros(11)

# 设置初始学习率为1，后续可以修改
learningRate = 1

2.4.4 寻找最优参数

现在我们可以使用和第一部分相同的优化函数来计算优化后的结果。

result2 = opt.fmin_tnc(func=costReg, x0=theta2, fprime=gradientReg, args=(X2, y2, learningRate))

得出的结果：(array([ 0.53010246, 0.29075567, -1.60725764, -0.58213819, 0.01781027,
-0.21329507, -0.40024142, -1.3714414 , 0.02264304, -0.95033581,
0.0344085 ]), 22, 1)

2.5 评估正则化逻辑回归

最后，我们可以使用第1部分中的预测函数来查看我们的方案在训练数据上的准确度。

theta_min = np.matrix(result2[0])
predictions = predict(theta_min, X2)
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y2)]
accuracy = (sum(map(int, correct)) % len(correct))
print ('accuracy = {0}%'.format(accuracy))

结果：accuracy = 78%

3、总结：

本次的试验还是挺有难度的，比如说这些理论和代码也不是一下就可以写出来的，想要写出这些代码需要的是对理论的烂熟于心和对python的熟悉。总的来说，这次的试验还是有很多没有看懂的地方，很多地方都是在依葫芦画瓢，之后还需要更多的时间来理解和消化。

4、试验用文件

百度网盘 提取码：1234