[Machine Learning] 线性回归与逻辑回归（含代码）

1 线性回归与逻辑回归的区别与联系

线性回归和逻辑回归都是广义上的线性回归。但两者又有不同。

1.1 决策边界

线性回归 是利用数理统计中的回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。表达形式为y = θx + e，其中e为误差服从均值为0的正态分布。适用于有监督学习的预测。

• 单变量线性回归：hθ(x)=θ0+θ1x，只包含一个自变量x和一个因变量hθ(x)，且二者的关系可以用一条直线近似表示；
• 多变量线性回归：hθ(x)=θ0+θ1x1+…+θnxn，包括两个或两个以上的自变量(x1、x2 …)，并且因变量和自变量是线性关系。

逻辑回归 首先通过sigmoid函数将样本映射到[0,1]之间的数值。sigmoid函数可以把任何连续的值映射到[0,1]之间，
决策边界：对多变量线性回归方程求sigmoid函数，

线性回归在整个实数域范围内进行预测，敏感度一致，而分类范围需要在[0,1]。逻辑回归是在线性回归的基础上，多一步sigmoid非线性映射。即减小预测范围，将预测值限定为[0,1]间的一种回归模型。所以线性回归模型对异常值敏感，而逻辑回归通过非线性变换减弱分离平面较远的点的影响，鲁棒性比线性回归的要好。

1.2 代价函数

代价函数是将随机事件或其有关随机变量的取值映射为非负实数以表示该随机事件的“风险”或“损失”的函数。用来评价模型的预测值和真实值不一样的程度，通常代价函数越小，模型的性能越好。

线性回归

• 单变量线性回归：
  

• 多变量线性回归：
  
  其中，

• m：训练样本的个数；

• hθ(x)：用参数θ和x预测出来的y值；

• y：原训练样本中的y值，也就是标准答案；

• 上角标(i)：第i个样本。

逻辑回归
对于线性回归，使用MSE得到的J(θ)为凸函数，但是对于logistic回归，由于进行了sigmoid非线性映射，在寻优时容易陷入局部最优的问题，所以考虑把sigmoid作log，得到的J(θ)。


J(θ)对θ求二阶导，得到其二阶导大于0，说明J(θ)为凸函数，使用梯度下降法寻优时，可以保证找到全局最小。

1.3 梯度下降法

使用梯度下降算法来求得能使代价函数最小的参数。

线性回归

逻辑回归

经梯度下降算法后线性回归与逻辑回归的公式形式看上去是一致的，但实际上两者是完全不同的，因为假设函数是不同的。
注：线性回归可使用正规方程求使代价函数最小的θ，但逻辑回归无法使用正规方程。

2 线性回归的代码实现

2.1 单变量线性回归

import numpy as np
import pandas as pd

# 读取数据
data = pd.read_csv(path, header=None, names=['x', 'y'])

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

# 处理数据及初始化
data.insert(0, 'Ones', 1) #在训练集中添加一列，以便可以使用向量化的解决方案来计算代价和梯度
cols = data.shape[1] #获取矩阵的列数
X = data.iloc[:,0:cols-1] #X是所有行，去掉最后一列
y = data.iloc[:,cols-1:cols] #y是所有行，最后一列

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

computeCost(X, y, 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

alpha = 0.01 #设置学习率α
iters = 2000 #设置迭代次数

g, cost = gradientDescent(X, y, theta, alpha, iters) # 运行梯度下降算法来将我们的参数θ适合于训练集
computeCost(X, y, g) #使用我们拟合的参数计算训练模型的代价函数

可视化结果

# 绘制线性模型
import matplotlib.pyplot as plt

x = np.linspace(data.x.min(), data.x.max(), 100)
f = g[0, 0] + (g[0, 1] * x)

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x, f, 'r', label='Prediction')
ax.scatter(data.x, data.y, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('xxx')
ax.set_ylabel('xxx')
ax.set_title('xxx')
plt.show()

# 绘制每个迭代次数与代价cost的关系图
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

2.2 多变量线性回归

import pandas as pd

data = pd.read_csv(path, header=None, names=['x1', 'x2', 'y'])

# 特征归一化
data = (data - data.mean()) / data.std()

# 第一列添加ones
data.insert(0, 'Ones', 1)

cols = data.shape[1]
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1]

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

# 调用梯度下降算法
g, cost = gradientDescent(X, y, theta, alpha, iters)

computeCost(X, y, g2)

除了调用梯度下降算法外，也可使用正规方程得到theta。

# 正规方程
def normalEqn(X, y):
    theta = np.linalg.inv(X.T@X)@X.T@y #X.T@X等价于X.T.dot(X)
    return theta

梯度下降与正规方程的比较

• 梯度下降：需要选择学习率α，需要多次迭代，当特征数量n大时也能较好适用，适用于各种类型的模型
• 正规方程：不需要选择学习率α，一次计算得出，但需要计算${\left(X^{T}X\right)}^{-1}$。如果特征数量n较大则运算代价大，因为矩阵逆的计算时间复杂度为$O(n3)$，通常来说当$n$小于10000 时还是可以接受的，只适用于线性模型，不适合逻辑回归模型等其他模型

2.3 sklearn

from sklearn import linear_model

model = linear_model.LinearRegression()
model.fit(X, y)

x = np.array(X[:, 1].A1)
f = model.predict(X).flatten()


fig, ax = plt.subplots(figsize=(12,8))
ax.plot( x, f, 'r', label='Prediction')
ax.scatter(data.x, data.y, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('xxx')
plt.show()

最小二乘法线性回归：sklearn.linear_model.LinearRegression(fit_intercept=True, normalize=False,copy_X=True, n_jobs=1)

主要参数说明

属性

• coef_：回归系数(斜率)

regr = linear_model.LinearRegression()
regr.fit(X,y)
regr.coef_

• intercept_：截距项

regr = linear_model.LinearRegression()
regr.fit(X,y)
regr.intercept_

主要方法

• fit(X, y, sample_weight=None)，x和y以矩阵的形式传入，sample_weight则是每条测试数据的权重，同样以矩阵方式传入
• predict(X)，预测方法，用来返回预测值
• score(X, y, sample_weight=None)，评分函数，将返回一个小于1的得分，可能会小于0
• get_params(deep=True)， 返回对regressor 的设置值

3 逻辑回归的代码实现

3.1 逻辑回归

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

data = pd.read_csv(path, header=None, names=['x1', 'x2', 'y'])

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

# 创建两个"正负"分数的散点图，并使用颜色编码来可视化
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['x1'], positive['x2'], s=50, c='b', marker='o', label='Yes')
ax.scatter(negative['x1'], negative['x2'], s=50, c='r', marker='x', label='No')
ax.legend()
ax.set_xlabel('x1 Score')
ax.set_ylabel('x2 Score')
plt.show()

# 定义sigmoid函数
def sigmoid(z):
    return 1 / (1 + np.exp(-z))

# 定义代价函数
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))

# 数据处理及初始化
data.insert(0, 'Ones', 1) #插入ones
cols = data.shape[1] #获得列数
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1:cols]

X = np.array(X.values)
y = np.array(y.values)
theta = np.zeros(3)

cost(theta, X, y)

# 梯度下降算法
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

gradient(theta, X, y) #查看数据和初始参数为0的梯度下降法的结果

# 用SciPy's truncated newton（TNC）实现寻找最优参数
import scipy.optimize as opt
result = opt.fmin_tnc(func=cost, x0=theta, fprime=gradient, args=(X, y))
cost(result[0], X, y)

# 使用获得的参数theta来为数据集X输出预测，并使用这个函数为分类器的训练精度打分。
def predict(theta, X):
    probability = sigmoid(X * theta.T)
    return [1 if x >= 0.5 else 0 for x in probability]

theta_min = np.matrix(result[0])
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))

若原始数据无明显的线性决策界限，一种方法是利用原始特征多项式得到特征，并对该特征使用像逻辑回归这样的线性技术。

degree = 5
x1 = data['Test 1']
x2 = data['Test 2']

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

for i in range(1, degree):
    for j in range(0, i):
        data['F' + str(i) + str(j)] = np.power(x1, i-j) * np.power(x2, j)

data.drop('Test 1', axis=1, inplace=True)
data.drop('Test 2', axis=1, inplace=True)

# 正则化代价函数
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

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

# 数据初始化，X、y
cols = data.shape[1]
X = data.iloc[:,1:cols]
y = data.iloc[:,0:1]

X = np.array(X.values)
y = np.array(y.values)
theta = np.zeros(11)

learningRate = 1 # 初始化学习率
costReg(theta, X, y, learningRate) #代价
gradientReg(theta, X, y, learningRate) #获得参数

result = opt.fmin_tnc(func=costReg, x0=theta, fprime=gradientReg, args=(X, y, learningRate))
theta_min = np.matrix(result[0])
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))

3.2 sklearn

from sklearn import linear_model #调用sklearn的线性回归包
model = linear_model.LogisticRegression(penalty='l2', C=1.0)
model.fit(X, y.ravel())

model.score(X, y)

参考资料

logistic 回归（内附推导）
线性回归和逻辑回归
浅析机器学习：线性回归 & 逻辑回归
代码部分主要引用黄海广博士的资料