2018-07-26-回归问题

2018-07-26 | 阅读：次

回归问题

假定我们现有一大批数据，包含房屋的面积和对应面积的房价信息，如果我们能得到房屋面积与房屋价格间的关系，那么，给定一个房屋时，我们只要知道其面积，就能大致推测出其价格了。

上面的问题还可以被描述为：

“OK，我具备了很多关于房屋面积及其对应售价的知识（数据），再通过一定的学习，当面对新的房屋面积时，我不再对其定价感到束手无策”。

通常，这类预测问题可以用回归模型（regression）进行解决，回归模型定义了输入与输出的关系，输入即现有知识，而输出则为预测。

一个预测问题在回归模型下的解决步骤为：

积累知识：我们将储备的知识称之为训练集 Training Set，很好理解，知识能够训练人进步
学习：学习如何预测，得到输入与输出的关系。在学习阶段，应当有合适的指导方针，江山不能仅凭热血就攻下。在这里，合适的指导方针我们称之为学习算法 Learning Algorithm
预测：学习完成后，当接受了新的数据（输入）后，我们就能通过学习阶段获得的对应关系来预测输出。

学习过程往往是艰苦的，“人谁无过，过而能改，善莫大焉”，因此对我们有这两点要求：

有手段能评估我们的学习正确性。
当学习效果不佳时，有手段能纠正我们的学习策略。

线性回归

预测

首先，我们明确几个常用的数学符号：

特征（feature）
特征向量（输入）
输出向量
假设（hypothesis）：

上面的表达式也称之为回归方程（regression equation）

误差评估

之前我们说到，需要某个手段来评估我们的学习效果。最常见的，我们通过最小均方（Least Mean Square）来描述误差：

误差评估的函数在机器学习中也称为代价函数（cost function）。

批量梯度下降

在引入了代价函数后，解决了“有手段评估学习的正确性”的问题，下面我们开始解决“当学习效果不佳时，有手段能纠正我们的学习策略”的问题。

在线性回归中，通常使用梯度下降（Gradient Descent）来调节 θ：

数学上，梯度方向是函数值下降最为剧烈的方向。那么，沿着 J(θ)J(θ) 的梯度方向走，我们就能接近其最小值，或者极小值，从而接近更高的预测精度。学习率 α是个相当玄乎的参数，其标识了沿梯度方向行进的速率，步子大了容易扯着蛋，很可能这一步就迈过了最小值。而步子小了，又会减缓我们找到最小值的速率。在实际编程中，学习率可以以 3 倍，10 倍这样进行取值尝试，如：

α=0.001,0.003,0.01…0.3,1

对于一个样本容量为 mm 的训练集，我们定义 θ的调优过程为：

我们称该过程为基于最小均方（LMS）的批量梯度下降法（Batch Gradient Descent），一方面，该方法虽然可以收敛到最小值，但是每调节一个 θjθj，都不得不遍历一遍样本集，如果样本的体积 mm很大，这样做无疑开销巨大。但另一方面，因为其可化解为向量型表示，所以就能利用到并行计算优化性能。

随机梯度下降

鉴于批量梯度下降的性能问题，又引入了随机梯度下降（Stochastic Gradient Descent）：

重复直到收敛（Repeat until convergence）

手段	概括	优点	缺点
批量梯度下降法	尽可能减小训练样本的总的预测代价	能够获得最优解，支持并行计算	样本容量较大时，性能显著下降
随机梯度下降法	尽可能的减小每个训练样本的预测代价	训练速度快	并不一定能获得全局最优，经常出现抖动和噪音，且不能通过并行计算优化

# coding: utf-8
# linear_regression/regression.py
import numpy as np
import matplotlib as plt
import time

def exeTime(func):
    """ 耗时计算装饰器
    """
    def newFunc(*args, **args2):
        t0 = time.time()
        back = func(*args, **args2)
        return back, time.time() - t0
    return newFunc

def loadDataSet(filename):
    """ 读取数据

    从文件中获取数据，在《机器学习实战中》，数据格式如下
    "feature1 TAB feature2 TAB feature3 TAB label"

    Args:
        filename 文件名

    Returns:
        X 训练样本集矩阵
        y 标签集矩阵
    """

    #统计特征数量
    numFeat = len(open(filename).readline().split('\t')) - 1
    X = []
    y = []
    #打开文件
    file = open(filename)
    for line in file.readlines():
        lineArr = []
        curLine = line.strip().split('\t')
        for i in range(numFeat):
            lineArr.append(float(curLine[i]))
        X.append(lineArr)
        y.append(float(curLine[-1]))
    return np.mat(X), np.mat(y).T

def h(theta, x):
    """预测函数
    
    Args:
        theta 相关系数矩阵
        x 特征向量

    Returns:
        预测结果
    """
    return (theta.T*x)[0,0]

def J(theta, X, y):
    """代价函数

    Args:
        theta 相关系数矩阵
        X 样本集矩阵
        y 标签集矩阵

    Returns:
        预测误差（代价）
    """
    m = len(X)
    return (X*theta-y).T*(X*theta-y)/(2*m)

@exeTime
def bgd(rate, maxLoop, epsilon, X, y):
    """批量梯度下降法

    Args:
        rate 学习率
        maxLoop 最大迭代次数
        epsilon 收敛精度
        X 样本矩阵
        y 标签矩阵

    Returns:
        (theta, errors, thetas), timeConsumed
    """
    m,n = X.shape
    # 初始化theta，n*1的矩阵
    theta = np.zeros((n,1))
    count = 0
    converged = False
    error = float('inf')
    errors = []
    thetas = {}
    for j in range(n):
        thetas[j] = [theta[j, 0]]
    while count <= maxLoop:
        if(converged):
            break
        count = count + 1
        for j in range(n):
            deriv = (y-X*theta).T*X[:, j]/m
            theta[j,0] = theta[j,0]+rate*deriv
            thetas[j].append(theta[j,0])
        error = J(theta, X, y)
        errors.append(error[0,0])
        # 如果已经收敛
        if(error < epsilon):
            converged = True
    return theta,errors,thetas

@exeTime
def sgd(rate, maxLoop, epsilon, X, y):
    """随机梯度下降法
    Args:
        rate 学习率
        maxLoop 最大迭代次数
        epsilon 收敛精度
        X 样本矩阵
        y 标签矩阵
    Returns:
        (theta, error, thetas), timeConsumed
    """
    m,n = X.shape
    # 初始化theta
    theta = np.zeros((n,1))
    count = 0
    converged = False
    error = float('inf')
    errors = []
    thetas = {}
    for j in range(n):
        thetas[j] = [theta[j, 0]]
    while count <= maxLoop:
        if converged:
            break
        count = count + 1
        errors.append(float('inf'))
        for i in range(m):
            if converged:
                break
            diff = y[i, 0]-h(theta, X[i].T)
            for j in range(n):
                theta[j,0] = theta[j,0] + rate*diff*X[i, j]
                thetas[j].append(theta[j, 0])
            error = J(theta, X, y)
            errors[-1] = error[0, 0]
            # 如果已经收敛
            if(error < epsilon):
                converged = True
    return theta, errors, thetas

def JLwr(theta, X, y, x, c):
    """局部加权线性回归的代价函数计算式

    Args:
        theta 相关系数矩阵
        X 样本集矩阵
        y 标签集矩阵
        x 待预测输入
        c tau
    Returns:
        预测代价
    """
    m,n = X.shape
    summerize = 0
    for i in range(m):
        diff = (X[i]-x)*(X[i]-x).T
        w = np.exp(-diff/(2*c*c))
        predictDiff = np.power(y[i] - X[i]*theta,2)
        summerize = summerize + w*predictDiff
    return summerize

@exeTime
def lwr(rate, maxLoop, epsilon, X, y, x, c=1):
    """局部加权线性回归

    Args:
        rate 学习率
        maxLoop 最大迭代次数
        epsilon 预测精度
        X 输入样本
        y 标签向量
        x 待预测向量
        c tau
    """
    m,n = X.shape
    # 初始化theta
    theta = np.zeros((n,1))
    count = 0
    converged = False
    error = float('inf')
    errors = []
    thetas = {}
    for j in range(n):
        thetas[j] = [theta[j,0]]
    # 执行批量梯度下降
    while count<=maxLoop:
        if(converged):
            break
        count = count + 1
        for j in range(n):
            deriv = (y-X*theta).T*X[:, j]/m
            theta[j,0] = theta[j,0]+rate*deriv
            thetas[j].append(theta[j,0])
        error = JLwr(theta, X, y, x, c)
        errors.append(error[0,0])
        # 如果已经收敛
        if(error < epsilon):
            converged = True
    return theta,errors,thetas

def standarize(X):
    """特征标准化处理

    Args:
        X 样本集
    Returns:
        标准后的样本集
    """
    m, n = X.shape
    # 归一化每一个特征
    for j in range(n):
        features = X[:,j]
        meanVal = features.mean(axis=0)
        std = features.std(axis=0)
        if std != 0:
            X[:, j] = (features-meanVal)/std
        else:
            X[:, j] = 0
    return X

def normalize(X):
    """特征归一化处理

    Args:
        X 样本集
    Returns:
        归一化后的样本集
    """
    m, n = X.shape
    # 归一化每一个特征
    for j in range(n):
        features = X[:,j]
        minVal = features.min(axis=0)
        maxVal = features.max(axis=0)
        diff = maxVal - minVal
        if diff != 0:
           X[:,j] = (features-minVal)/diff
        else:
           X[:,j] = 0
    return X

网上例子：

https://blog.csdn.net/SA14023053/article/details/51703204披萨的直径与价格

#导入线性模型
from sklearn import linear_model
#绘图
import matplotlib.pyplot as plt
import numpy as np


#X表示匹萨尺寸 Y表示匹萨价格
X = [[6], [8], [10], [14], [18]]
Y = [[7], [9], [13], [17.5], [18]]

print(X)
print(Y)

#使用线性回归
clf = linear_model.LinearRegression()
#导入数据集fit(x, y): 训练。分析模型参数，填充数据集。其中x为特征，y位标记或类属性。
#predict(): 预测。它通过fit()算出的模型参数构成的模型，对解释变量进行预测其类属性。预测方法将返回预测值y_pred。
clf.fit(X, Y)

#预测数据集
print(np.array([12]).reshape(-1, 1))
res = clf.predict(np.array([12]).reshape(-1, 1))[0]
print(res)

#预测结果
X2 = [[0], [10], [14], [25]]
Y2 = clf.predict(X2)
print(Y2)

#预测结果
X3 = [[0]]
Y3 = clf.predict(X3)
print(Y3)

#绘图
plt.figure()
plt.title('1')
plt.xlabel('x')
plt.ylabel('cost')
plt.axis([0, 25, 0, 25])
plt.grid(True)
#绘制训练集散点图
plt.plot(X, Y, 'k.')
#绘制预测集直线
plt.plot(X2, Y2, 'g-')
plt.show()

运行结果：

D:\1b\Anoconda\setup\set\python.exe E:/python_workspace/ML/WuEnDa_Work/LInerReg.py
[[6], [8], [10], [14], [18]]
[[7], [9], [13], [17.5], [18]]
[[12]]
[ 13.68103467]
[[  1.96551743]
 [ 11.72844846]
 [ 15.63362088]
 [ 26.37284501]]
[[ 1.96551743]]

标准化实验：

import numpy as np

def standardize(X):
    #m表示行，n表示列
    m ,n = X.shape
    for j in range(n):
        features = X[:, j]
        meanVal = features.mean(axis = 0)
        std = features.std(axis = 0)
        print("std", std)
        print("meanVal", meanVal)
        print("feature", features)
        if std != 0:
            X[:, j] = (features - meanVal)/std
        else:
            X[:, j] = 0
    return X

if __name__ == "__main__":
    X = np.arange(0, 9).reshape(9, 1)
    #X = [[0], [10], [14], [25]]
    s = standardize(X)
    print("S", s)