局部加权线性回归和岭回归

"""
线性回归的一个问题就是有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏
估计，显而易见，如果模型欠拟合将不能取得最好的预测效果，所以有些方法允许在估计中引入一些
偏差，从而降低预测的均方误差
"""
from numpy import *
import matplotlib.pyplot as plt


def loadDataSet(fileName):
    """
    该函数打开一个用tab键分隔的文本文件，这里仍然默认文件每行的最后一个值是目标值
    readline():该方法每次读取一行内容，所以读取时占用内存小，比较适合大文件，该方法返回一个字符串对象
    readlines():该方法读取整个文件的所有行，保存在一个列表变量中，每行作为一个元素，但读取大文件比较占用内存
    该函数第一个返回值返回列表中的自变量第一列和第二列，　第二个返回值返回结果，最后一列
    :param fileName:
    :return:
    """
    numFeat = len(open(fileName).readline().split("\t")) - 1
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():  # readlines()读取整个文件保存在列表变量中
        lineArr = []
        curLine = line.strip().split("\t")  # strip()返回移除字符串头尾制定的字符生成的新
        for i in range(numFeat):
            lineArr.append(float(curLine[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    return dataMat, labelMat


def standRegres(xArr, yArr):
    """
    该函数用来计算最佳拟合直线，函数首先读入x和y，并将他们保存到矩阵中，然后计算x的转置与x的乘积
    然后判断他的行列式是否为零，如果行列式为零，那么计算逆矩阵就会出现错误NumPy提供一个线性代数
    的库linalg，可以通过linalg.det()来计算行列式，最后，如果行列式非零，计算并返回w
    :param xArr:
    :param yArr:
    :return:
    """
    xMat = mat(xArr)  # 列表类型转换为矩阵类型
    yMat = mat(yArr).T
    xTx = xMat.T * xMat  # 求x的转置与x相乘的结果

    # 判断矩阵的行列式是否为零，如果为零，则无法计算逆矩阵
    if linalg.det(xTx) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T * yMat)  # xTx.I 为xTx的逆，计算得出未知量的系数
    return ws


def lwlr(testPoint, xArr, yArr, k = 1.0):
    """
    该函数的功能是给定x空间中的任意一点，计算出对应的预测值yHat,函数lwlr()的开头
    读入数据并创建所需要的矩阵，之后创建对角权重矩阵weights。权重矩阵是一个方阵，阶数等于
    样本点个数，也就是说，该矩阵为每个样本点初始化了一个权重，接着遍历数据集，
    计算每个样本点的权重值，随着样本点与待预测点距离的递增，权重将以指数级递减
    输入参数k控制衰减的速度。在权重矩阵计算完毕后，就可以得到对回归系数ws的一个估计
    :param testPoint:
    :param xArr:
    :param yArr:
    :param k:
    :return:
    """
    xMat = mat(xArr) # 列表类型转换为矩阵类型
    yMat = mat(yArr).T
    m = shape(xMat)[0] # shape()功能是读取矩阵的长度
    weights = mat(eye(m)) # 生成m阶单位阵
    for j in range(m):
        diffMat = testPoint - xMat[j, :]
        weights[j, j] = exp(diffMat * diffMat.T / (-2.0 * k ** 2))
    xTx = xMat.T * (weights * xMat)
    if linalg.det(xTx) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T * (weights * yMat))
    return testPoint * ws


def lwlrTest(testArr, xArr, yArr, k=1.0):
    m = shape(testArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i], xArr, yArr, k)
    return yHat


def rssError(yArr, yHatArr):
    """
    计算误差
    :param yArr:
    :param yHatArr:
    :return:
    """
    return ((yArr-yHatArr)**2).sum() # (实际值－预测值)**2再求和

"""
岭回归：最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计
这里通过引入　来限制所有ｗ之和，通过引入该惩罚项能够减少不重要的参数这个技术在统计学
中也叫做缩减
缩减法可以去掉不重要的参数，因此能更好的理解数据，此外，与简单的线性回归相比
缩减法能取得更好的预测效果
""",


def ridgeRegres(xMat, yMat, lam=0.2):
    """
    该函数用于计算回归系数，函数实现了给定lambda下的岭回归求解，如果没有指定lambda
    则默认为0.2,由于lambda是Python保留的关键字，因此程序使用了lam来代替，该函数首先
    构建矩阵xTx　然后用lam乘以单位矩阵(使用NumPy库中的eye方法生成)在普通回归方法可
    能会产生错误时，岭回归仍然可以正常工作
    :param xMat:
    :param yMat:
    :param lam:
    :return:
    """
    xTx = xMat.T * xMat
    denom = xTx + eye(shape(xMat)[1]) * lam
    if linalg.det(denom) == 0:
        print("This matrix is singular, cannot so inverse")
        return
    ws = denom.I * (xMat.T * yMat)
    return ws

def draw():
    abX, abY = loadDataSet("abalone.txt")
    yHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)
    yHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)
    yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)
    r1 = rssError(abY[0:99], yHat01.T)
    r2 = rssError(abY[0:99], yHat1.T)
    r3 = rssError(abY[0:99], yHat10.T)
    print(r1)
    print(r2)
    print(r3)
    yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
    rr1 = rssError(abY[100:199], yHat01.T)
    yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)
    rr2 = rssError(abY[100:199], yHat1.T)
    yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)
    rr3 = rssError(abY[100:199], yHat10.T)
    print(rr1)
    print(rr2)
    print(rr3)


def ridgeTest(xArr, yArr):
    """
    该函数用于在一组lambda上测试结果
    为了使用岭回归和缩减技术，首先需要对特征做标准化处理，
    使每维特征具有相同的重要性(不考虑特征代表什么)
    该函数展示了标准化的过程，具体的做法是所有特征都减去各自的均值并除以方差
    :param xArr:
    :param yArr:
    :return:
    """
    xMat = mat(xArr) # 转换为矩阵
    yMat = mat(yArr).T # 转换为矩阵，并进行转置

    # 进行标准化操作
    yMean = mean(yMat, 0) # 对所有元素求和再求出均值
    yMat = yMat - yMean
    xMeans = mean(xMat, 0)
    xVar = var(xMat, 0) # 求xMat的方差
    xMat = (xMat - xMeans) / xVar # 进行标准化处理
    numTestPts = 30
    wMat = zeros((numTestPts, shape(xMat)[1])) # 生成行数为numTestPts,列数为shape(xMat)[1]的零矩阵
    for i in range(numTestPts):
        ws = ridgeRegres(xMat, yMat, exp(i - 10))
        # print(ws)
        wMat[i, :] = ws.T
    return wMat


def regularize(xMat):
    inMat = xMat.copy()
    inMeans = mean(inMat, 0)
    inVar = var(inMat, 0)
    inMat = (inMat - inMeans) / inVar
    return inMat


def stageWise(xArr, yArr, eps = 0.01, numIt = 100):
    """
    该函数是一个逐步线性回归算法的实现，它与lasso做法相近，但是计算简单
    该函数的输入包括：输入数据xArr和预测变量yArr。此外还有两个参数，一个
    是esp,表示每次迭代需要调整的步长，numIt表示需要调整的次数
    :param xArr:
    :param yArr:
    :param eps:
    :param numIt:
    :return:
    """
    # 输入数据转换并存入矩阵中，然后把特征按照均值为０方差为１进行标准化处理
    xMat = mat(xArr)
    yMat = mat(yArr).T
    yMean = mean(yMat, 0)
    yMat = yMat - yMean
    xMat = regularize(xMat)
    m, n = shape(xMat) # 返回xMat矩阵的行数和列数
    returnMat = zeros((numIt, n))

    # 创建向量ws来保存w的值，并且为了实现贪心算法建立了ws的两份副本，接下来
    # 优化过程中需要迭代numIt次，并且在每次迭代时都打印出w向量，用于分析算法执行过程和效果
    ws = zeros((n, 1)) # 生成n行1列的零矩阵
    wsTest = ws.copy()
    wsMax = ws.copy()
    for i in range(numIt):
        print(ws.T)
        lowestError = inf # 误差初始值设置为正无穷
        for j in range(n):
            for sign in [-1, 1]:
                wsTest = ws.copy()
                wsTest[j] += eps * sign
                yTest = xMat * wsTest
                rssE = rssError(yMat.A, yTest.A)
                if rssE < lowestError:
                    lowestError = rssE
                    wsMax = wsTest
        ws = wsMax.copy()
        returnMat[i, :] = ws.T
    return returnMat


if __name__ == '__main__':
    abX, abY = loadDataSet("abalone.txt")
    stageWise(abX, abY)
    ridgeWeights = ridgeTest(abX, abY)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(ridgeWeights)
    plt.show()