前言

​ 线性回归只能拟合线性曲面（广义的曲面），如果一个回归任务中的输出变量 $y\left(y\in \mathbb{R}\right)$ 关于特征向量 $x=\left(x_0,x_1,\cdots ,x_n\right)\left(x\in {\mathbb{R}}^{n+1},x_0=1\right)$ 不是线性的，那怎样拟合可以比较好呢？局部加权线性回归是一种方法。

1. 局部加权线性回归模型

​ 当两个样本点的特征向量距离较近时，它们的输出变量通常比较接近，也就是说它们具有相似的性质。现在我们有 $m$ 个训练样本点 $\left(x^{\left(1\right)},y^{\left(1\right)}\right),\left(x^{\left(2\right)},y^{\left(2\right)}\right),\cdots ,\left(x^{\left(m\right)},y^{\left(m\right)}\right)$ 构成训练集。那么当我们给一个特征向量 $x$ 时，要根据训练集中的信息预测出它的输出变量 $y$ ，距离 $x$ 越近的训练集中的样本点越应该被重视。换句话说，我们需要给训练集中的每个样本点一个权值，距离 $x$ 越近的样本点的权值应该越大，这是赋权的一个原则，那么具体应该如何赋权呢？

​ 首先应该想一想，我们怎样衡量距离呢？欧氏距离（2范数）是一种常用的衡量方式，即：
$$dist\left(x,x^{\left(i\right)}\right)={\Vert x-x^{\left(i\right)}\Vert }_2$$
我们用下面的公式给训练集中第 $i$ 样本点赋权值：
$$w_i=\exp \left(-\frac{{\Vert x-x^{\left(i\right)}\Vert }_2^2}{2{\sigma }^2}\right)$$
这里似乎很突兀，还不知道为什么，直接就给出了这样的公式，我也不太清楚为什么，这个算法经常使用这样的函数来赋权，对于不是很想深入理论研究的我姑且认为是玄学吧！当然肯定不是玄学，这样做应该是有一定理论依据的，不过我就暂时不了解了。

​

​ 对于给定的特征向量作为输入，我们要预测它的输出变量，使用下面的经过加权处理后的损失函数：
$$J\left(\theta \right)=\frac{1}{2}\sum _{i=1}^m{w}_i{\left(h_{\theta }\left(x^{\left(i\right)}\right)-y^{\left(i\right)}\right)}^2$$
这里还是挺自然的，就是一个简单的加权求和，最小化这样的损失函数时，算法肯定会寻找使得距离 $x$ 越近的样本点的损失越小的 $\theta$ ，也就是说算法更加重视了那些距离 $x$ 近的点，这也有利于更准确地预测 $y$ 。该函数使用矩阵的方式写出来就是：
$$J\left(\theta \right)=\frac{1}{2}{\left(X\theta -y\right)}^{T}W\left(X\theta -y\right)$$
其中权值矩阵 $W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(w_1,w_2,\cdots ,w_m\right)$ 。

​ 优化这个损失函数的过程就是该算法的训练过程。优化这个损失函数后，就可以使用得到参数 $\theta$ ，进而使用 $h_{\theta }\left(x\right)$ 预测输出变量 $y$ 。 如果我们有 $K$ 个特征向量需要被预测，记这 $K$ 个特征向量为 $x^{\left(m+1\right)},x^{\left(m+2\right)},\cdots ,x^{\left(m+K\right)}$ ，那么对于每一个特征向量 $x^{\left(m+i\right)}\left(1⩽i⩽K\right)$ 我们都要计算一次权值矩阵，然后然后优化损失函数得到 $\theta$ ，进而预测输出变量 $y^{\left(m+i\right)}$ 。

2. 求解方法

​ 梯度下降法是一种很常用的优化方法，这里当然也可以使用梯度下降法，主要就是求解出损失函数的梯度，求解梯度实际上也就是求导数，对其求导：
$$\frac{\partial J\left(\theta \right)}{\partial \theta }=X^{T}W\left(X\theta -y\right)=X^{T}WX\theta -X^{T}Wy$$
这样也就可以使用梯度下降法了，这里就不多说了，直接看正规方程的方法，令其导数为 $0$ 可以解出该函数的极小值点：
$$\theta ={\left(X^{T}WX\right)}^{-1}{X}^{T}Wy$$
​ 值得注意的是，局部加权线性回归时间代价十分高，因为对于每一个需要预测的点，我们都要求出相对应的权值矩阵并且单独训练。而普通的线性回归算法只需要训练一次即可预测所有的点。

​ 局部加权线性回归的好处就是它几乎可以拟合任意形状，只要超参数 $\sigma$ 选取合适，那么该算法可以拟合的较好。实际生活中，数据的输出变量和特征向量是线性关系的情况不多，因此纯线性回归模型可能适用性不强，但是局部加权线性回归可以对非线性的数据分布做出较好的拟合。

3. 代码实现

​ 写出局部加权线性回归函数。

import numpy as np

def calweight(a, b): # 求解权值的函数
    sigma = 0.3
    return np.exp(-np.linalg.norm(a-b)**2 / (2*sigma**2))

def LWLR(X_train, y_train, X_test): # 局部加权线性回归
    X_train = np.insert(X_train,0,1,axis=1)
    X_test = np.insert(X_test, 0, 1, axis = 1)
    y_pred = np.empty(X_test.shape[0]) # 预测值
    
    for k in range(X_test.shape[0]): # 对所有的需要预测的特征向量先训练再预测
        # 求解权值矩阵   
        w = np.empty(X_train.shape[0])
        for i in range(X_train.shape[0]):
            w[i] = calweight(X_test[k], X_train[i])
        W = np.diag(w) #权值矩阵
        
        theta=np.dot(np.linalg.pinv(np.dot(np.dot(X_train.T,W),X_train)),np.dot(np.dot(X_train.T,W),y_train))
        y_pred[k] = np.dot(theta,X_test[k]) # 预测
        
    return y_pred

​ 创造数据，使用该函数拟合回归曲线。

import matplotlib.pyplot as plt

def CreateData():
    X = np.arange(0,10,0.1)
    y = np.empty(X.shape[0])
    for i in range(X.shape[0]):
        y[i] = X[i]**3 - 10*X[i]**2 + X[i] + np.random.uniform(-10,10)
    return X[:,np.newaxis], y

np.random.seed(0)
X, y = CreateData()
plt.scatter(X, y) # 可视化
X_test = np.arange(0.05,10,0.1) #测试集，目的是通过其画出拟合曲线
X_test = X_test[:,np.newaxis]
y_pred = LWLR(X, y, X_test)
plt.plot(X_test, y_pred, color='red') # 画出拟合曲线
plt.show()

​ 拟合曲线如下：

​ $\sigma =2$ ，欠拟合

​ $\sigma =1$ ，感觉在底部稍微有点欠拟合

​ $\sigma =0.5$ ，拟合效果不错

​ $\sigma =0.1$ ，有些过拟合

​ 只要超参数 $\sigma$ 选取合适，很多情况下都会有较好的拟合效果。

局部加权线性回归

具体理论见上次笔记《线性回归》

预测鲍鱼年龄

import numpy as np

class LocalWeightedLinearRegression(object):
    def __init__(self,train_data,train_result):
        """
        :param train_data: 输入训练数据
        :param train_result: 训练数据真实结果
        """
        # 获得输入数据集的形状
        row, col = np.shape(train_data)
        # 构造输入数据数组
        self.Train_Data = [0] * row
        # 给每组输入数据增添常数项1
        for (index, data) in enumerate(train_data):
            Data = [1.0]
            # 把每组data拓展到Data内，即把每组data的每一维数据依次添加到Data
            Data.extend(list(data))
            self.Train_Data[index] = Data
        self.Train_Data = np.array(self.Train_Data)
        # 构造输入数据对应的结果
        self.Train_Result = train_result
        # 定义数据权重
        self.weight = np.zeros((row,row))
        # 定义局部加权回归模型参数
        self.Theta = []

    def Gaussian_Weight(self,data,k):
        """
        这是计算测试权重的函数
        :param data: 输入数据
        :param k: 带宽系数
        """
        # data的数据类型是np.array，那么利用dot方法
        # 进行矩阵运算的结果是矩阵，哪怕只有一个元素
        sum = np.sum(data*data)
        return np.exp(sum/(-2.0*k**2))

    def predict_NormalEquation(self,test_data,k):
        """
        这是利用正规方程对测试数据集的局部加权线性回归预测函数
        :param test_data: 测试数据集
        :param k: 带宽系数
        """
        # 对测试数据集全加入一维1，以适应矩阵乘法
        data = []
        for test in test_data:
            # 对测试数据加入1维特征，以适应矩阵乘法
            tmp = [1.0]
            tmp.extend(test)
            data.append(tmp)
        test_data = np.array(data)
        # 计算test_data与训练数据集之间的权重矩阵
        for (index, train_data) in enumerate(self.Train_Data):
            diff = test_data-self.Train_Data[index]
            self.weight[index, index] = self.Gaussian_Weight(diff, k)
        # 计算XTWX
        XTWX = self.Train_Data.T.dot(self.weight).dot(self.Train_Data)
        """
          0.001*np.eye(np.shape(self.Train_Data.T))是
          防止出现原始XT的行列式为0，即防止原始XT不可逆
       """
        # 获得输入数据数组形状
        row, col = np.shape(self.Train_Data)
        # 若XTWX的行列式为0，即XTWX不可逆，对XTWX进行数学处理
        if np.linalg.det(XTWX) == 0.0:
            XTWX = XTWX + 0.001 * np.eye(col, col)
        # 计算矩阵的逆
        inv = np.linalg.inv(XTWX)
        # 计算模型参数Thetha
        XTWY = self.Train_Data.T.dot(self.weight).dot(np.reshape(self.Train_Result,(len(self.Train_Result),1)))
        self.Theta = inv.dot(XTWY)
        # 对测试数据test_data进行预测
        predict_result = test_data.dot(self.Theta).T[0]
        return predict_result

import LocalWeightedLinearRegression as LWLR
from sklearn.model_selection import train_test_split
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

def Load_Abalone(path):
    """
    这是导入鲍鱼数据集的函数
    :param path: 文件路径
    """
    # 定义鲍鱼数据和结果数组
    data = []
    result = []
    # 打开路径为path的文件
    with open(path) as f:
        # 遍历文件中的每一行
        for line in f.readlines():
            str = line.strip().split(',')
            tmp = []
            length = len(str[1:])
            # 鲍鱼数据集中的第一个属性是性别，该属性属于离散数据
            # 因此导入数据时必须抛开这一列，最后一列是环数，加1.5可以预测年龄
            for (index,s) in enumerate(str[1:]):
                # 最后一个数据追加到result
                if index == length-1:
                    result.append(float(s)+1.5)
                # 否则剩下的数据追加到tmp临时数组
                else:
                    tmp.append(float(s))
            # 一组数据追加到数据集中
            data.append(tmp)
        data = np.array(data)
        result = np.array(result)
    return data,result

def Merge(data,col):
    """
    这是生成DataFrame数据的函数
    :param data:输入数据
    :param col:列名称数组
    """
    Data = np.array(data).T
    return pd.DataFrame(Data,columns=col)

def run_main():
    """
       这是主函数
    """
    # 导入鲍鱼数据集
    path = "./abalone_data.txt"
    Data,Result = Load_Abalone(path)

    # 把数据集分成训练集合测试集
    Train_Data,Test_Data,Train_Result,Test_Result = train_test_split\
        (Data,Result,test_size=0.01,random_state=10)

    # 解决Matplotlib中的中文乱码问题，以便于后面实验结果可视化
    mpl.rcParams['font.sans-serif'] = [u'simHei']
    mpl.rcParams['axes.unicode_minus'] = False

    # 可视化测试集
    col = ['长度','直径','高度','总重量','剥壳重量','内脏重量','壳重']

    # 遍历局部加权线性回归算法的预测结果
    fig = plt.figure()
    for (index, c) in enumerate(col):
        if index == 0:
            ax = fig.add_subplot(311)
        else:
            ax = fig.add_subplot(334+index-1)
        ax.scatter(Test_Data[:,index], Test_Result, alpha=0.5, c='b', s=10)
        ax.grid(True)
        plt.xlabel(c)
        plt.ylabel("鲍鱼年龄")
        # 子图之间使用紧致布局
        plt.tight_layout()
    plt.savefig("./测试结果可视化.jpg", bbox_inches='tight')
    plt.close()

    for (index,c) in enumerate(col):
        plt.scatter(Test_Data[:,index],Test_Result,alpha=0.5,c='b',s=10)
        plt.grid(True)
        plt.xlabel(c)
        plt.ylabel("鲍鱼年龄")
        plt.savefig("./"+c+"可视化.jpg",bbox_inches='tight')
        #plt.show()
        plt.close()

    # 初始化局部加权线性回归模型
    lwlr = LWLR(Train_Data,Train_Result)

    # 初始化局部加权回归带宽系数参数
    K = [0.0001,0.001,0.003,0.005,0.01,0.05,0.1,0.3,0.5]
    tmp = list(np.arange(1,101))
    K.extend(tmp)
    predict_result = []
    Loss = []
    # 把测试集进行从小到大排序
    sort = np.argsort(Test_Data[:, 1])
    Test_Data = Test_Data[sort]
    Test_Result = Test_Result[sort]
    # 遍历每个带宽系数，利用局部加权线性回归算法进行预测,并计算预测误差
    for k in K:
        # 在带宽系数k下，利用局部加权线性回归进行预测
        predict = lwlr.predict_NormalEquation(Test_Data,k)
        #print(np.shape(predict))
        # 计算每组数据的预测误差
        loss = (Test_Result-predict)**2
        print("k=%f时的误差：%f"%(k,np.sum(loss)))
        predict_result.append(predict)
        Loss.append(loss)
        # 可视化预测结果
        plt.scatter(Test_Data[:, 1], Test_Result, alpha=0.5, c='b', s=10)
        plt.plot(Test_Data[:,1],predict,'r')
        plt.grid(True)
        plt.xlabel('直径')
        plt.ylabel("鲍鱼年龄")
        plt.savefig("./k=" + str(k) + "可视化.jpg", bbox_inches='tight')
        plt.close()

    # 部分预测结果可视化
    # k = [0.1,0.3,1,3,10,100]
    k = [0.1,0.3,1,3,10,100]
    index= [6,7,9,11,18,108]
    # 遍历每个带宽系数，利用局部加权线性回归算法进行预测,并计算预测误差
    fig = plt.figure()
    for (j,(i,k_)) in enumerate(zip(index,k)):
        # 在带宽系数k下，利用局部加权线性回归进行预测
        predict = predict_result[i]
        # 可视化预测结果
        ax = fig.add_subplot(230+j+1)
        ax.scatter(Test_Data[:, 1], Test_Result, alpha=0.5, c='b', s=10)
        ax.plot(Test_Data[:,1],predict,'r')
        ax.grid(True)
        ax.legend(labels=["k=" + str(k_)], loc="best")
        plt.xlabel('直径')
        plt.ylabel("鲍鱼年龄")
        plt.tight_layout()
    plt.savefig("./部分预测结果.jpg")

    # 保存预测数据
    data = [Test_Result]
    col = ["真实结果"]
    for (index,k) in enumerate(K):
        data.append(predict_result[index])
        col.append("k="+str(k))
    Data = Merge(data,col)
    Data.to_excel("./真实结果与预测结果.xlsx")

    # 保存预测误差结果及其统计信息
    data = []
    col = []
    for (index, k) in enumerate(K):
        data.append(Loss[index])
        col.append("k=" + str(k))
    Data = Merge(data,col)
    Data.to_excel("./预测误差.xlsx")
    information = Data.describe()
    information.to_excel("./预测误差统计信息.xlsx")

    # 可视化不同带宽参数的局部加权线性回归模型在测试集的均方误差和预测标准差
    K = list(np.arange(1, 101))
    col = ["LWLR-MSE", "LWLR-std"]
    LWLR_MSE = list(information.loc['mean'])[9:]
    LWLR_std = list(information.loc['std'])[9:]
    plt.plot(K, LWLR_MSE, 'b')
    plt.plot(K, LWLR_std, 'c-.')
    plt.grid(True)
    plt.legend(labels=col, loc='best')
    plt.xlabel("带宽系数")
    plt.savefig("./局部加权线性回归的预测均方误差和标准差.jpg", bbox_inches='tight')
    plt.show()
    plt.close()

if __name__ == '__main__':
    run_main()

线性回归虽然简单，但是容易出现问题：就是“欠拟合”和“过拟合”，欠拟合是由于我们并不能很好的拟合我们的训练数据，导致出现较大的训练误差；而过拟合是由于我们过度拟合训练数据，导致我们的模型过度复杂而产生较大的测试误差。

如下图所示：左边的图就是欠拟合，很明显我们用直线是无法很好的拟合训练数据的，最右边的就是产生了过拟合，中间的曲线就是比较好的。

解决欠拟合我们可以增加额外的特征，或者增加多项式（如$x^2$之类的）。

接下来介绍一种解决欠拟合的方法——局部加权线性回归。

1.代价函数：
$$J(\theta )=\sum _{i=1}^m{w}_{(i,i)}(f(x_i)-y_i{)}^2$$
其中$f(x_i)$是我们的预测值，$f(x_i)={\theta }^T{x}_i$；$w_{(i,i)}$是权重，它是通过要预测的点i与数据集中点的距离来确定的，距离越近， 权值越大，相应的误差值$(f(x_i)-y)$影响就越大（可以达到增加多项式的效果），反之同理。我们通常使用如下的权重函数：
$$w_{(i,i)}=exp(-\frac{(x_i-x{)}^2}{2k^2})$$
其中，x是我们的要预测的点，k是我们需要指定的参数，他控制了权值随距离变化的速率。
接下来我们可以通过梯度下降或者正规方程的方法来求解我们的$\theta$，这里使用正规方程的解法，首先将我们的上面的代价函数写成矩阵的形式：
$$J(\theta )=(X\theta -y{)}^{T}W(X\theta -y)$$
其中X是我们的训练集，每一行是一个样本，y是我们的真实值，W是m大小的方阵，并且只有对角线上有值，其他地方都为0，（这里不懂得可以随便去一个样本大小为例子乘以下就知道了）。

最后我们就用对$J(\theta )$关于$\theta$求导，并令其为0求出我们的$\theta$等于：
$$\theta =(X^{T}WX{)}^{-1}{X}^{T}Wy$$

2.代码实现：
这里直接用了机器学习实战上面的实现代码
2.1算法实现：

from numpy import *
def lwlr(testPoint,xArr,yArr,k=1.0):
    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 testPoint * ws

def lwlrTest(testArr,xArr,yArr,k=1.0):  #testArr是我们要预测的所有点的矩阵
    m = shape(testArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)   #对于每一个要预测的点都要重新计算权值
    return yHat

2.2绘制结果：

from numpy import *
import matplotlib.pyplot as plt
import regression

xArr, yArr = regression.loadDataSet('data/Ch08/ex0.txt')
yHat = regression.lwlrTest(xArr,xArr,yArr,1.0)
xMat = mat(xArr)
srtInd = xMat[:,1].argsort(0)
xSort = xMat[srtInd][:,0,:]

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xSort[:,1],yHat[srtInd])
ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0],s=2,c='red')
plt.show()

3.运行结果：
k=1.0（欠拟合）：

k=0.01（拟合的比较好）：

k=0.003（过拟合）：

在调用上面代码时，我们通过会选取不同的k值，以获得最佳效果。
通过上面我们也知道每预测一个点，我们都要重新计算权重W矩阵，增加了我们的计算量。