回归预测（7）

当我们拥有2021年每月降水量，想要获得2022年每月降水量时，首先想到的就是估计降水量与各月的函数关系，从而预测未知月份的降水量。

而解决上述问题，首先想到的算法就是--最小二乘

构建如y=Bx的函数关系，可以展开成y=b0*cosx+b1*x+b2*tanx+……,x为解决该问题的基函数，B为对应每项基函数的系数阵，而最小二乘求解的方式也很简单：

在求解得到系数阵后，对于新得到的x就能成功预测新的y值，代码实现如下：

def standRegres(xArr,yArr):
    #输入：xArr--基函数值，yArr--预测对象值
    xMat = mat(xArr); yMat = mat(yArr).T
    xTx = xMat.T*xMat
    if linalg.det(xTx) == 0.0:
        print("This matrix is singular, cannot do inverse")    #奇异矩阵
        return
    ws = xTx.I * (xMat.T*yMat)
    return ws

而简单最小二乘没法解决所有问题，例如在影响因子A,B,C,D中，可能待预测目标受A的影响最大，而D对目标的影响很小，采用简单最小二乘对其一概而论就不再适用。这时就需要引入权重因子P，对每个影响因子的干扰定不同的权重P(A)、P(B)、P(C)、P(D)，则函数模型变为：

 而权重可以通过不同的权函数来进行计算获取，代码实现如下：

def lwlr(xArr,yArr,k=1.0):
    #输入：xArr--影响因素；yArr--预测目标；k--权函数参数
    xMat = mat(xArr); yMat = mat(yArr).T
    m = shape(xMat)[0]
    weights = mat(eye((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 ws

但当影响因素的个数大于样本个数时，即该样本无法完成影响因素的描述，换言之系数阵秩亏。则需要去除几个与样本关系不紧密的影响因素，也就是缩减系数，通常的方法有岭回归。也就是在简单最小二乘中引入岭参数，计算公式如下：

代码实现如下：

def ridgeRegres(xMat,yMat,lam=0.2):
    #输入：xMat、yMat--与之前相似，lam--岭参数
    xTx = xMat.T*xMat
    denom = xTx + eye(shape(xMat)[1])*lam
    if linalg.det(denom) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = denom.I * (xMat.T*yMat)
    return ws

除了岭回归外，还可以采用--前向逐步回归

迭代确定权值，直至误差最小。代码实现如下：

def stageWise(xArr,yArr,eps=0.01,numIt=100):
    #输入：xArr,yArr--自变量与预测目标；eps--容许误差；numIt--迭代次数
    xMat = mat(xArr); yMat=mat(yArr).T
    yMean = mean(yMat,0)
    yMat = yMat - yMean
    xMat = regularize(xMat)    #重心化处理
    m,n=shape(xMat)
    returnMat = zeros((numIt,n))
    ws = zeros((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