局部加权回归r语言代码 局部加权线性回归

一、问题引入

我们现实生活中的很多数据不一定都能用线性模型描述。依然是房价问题，很明显直线非但不能很好的拟合所有数据点，而且误差非常大，但是一条类似二次函数的曲线却能拟合地很好。为了解决非线性模型建立线性模型的问题，我们预测一个点的值时，选择与这个点相近的点而不是所有的点做线性回归。基于这个思想，便产生了局部加权线性回归算法。在这个算法中，其他离一个点越近，权重越大，对回归系数的贡献就越多。

二、问题分析

本算法依然使用损失函数J，只不过是加权的J函数：

其中w(i)是权重，它根据要预测的点与数据集中的点的距离来为数据集中的点赋权值。当某点离要预测的点越远，其权重越小，否则越大。一个比较好的权重函数如下：

该函数称为指数衰减函数，其中k为波长参数，它控制了权值随距离下降的速率，该函数形式上类似高斯分布(正态分布)，但并没有任何高斯分布的意义。该算法解出回归系数如下：

在使用这个算法训练数据的时候，不仅需要学习线性回归的参数，还需要学习波长参数。这个算法的问题在于，对于每一个要预测的点，都要重新依据整个数据集计算一个线性回归模型出来，使得算法代价极高。

三、代码实现
1.Matlab版
(1) 局部加权线性回归函数

% 对所有点计算估计值
function yHat = lwlr(testMat, xMat, yMat, k)
    [row, ~] = size(testMat);
    yHat = zeros(1, row);
    for i = 1:1:row
        yHat(i) = lwlrPoint(testMat(i,:), xMat, yMat, k);
    end
end
% 对单个点计算估计值
function yHatPoint = lwlrPoint(point, xMat, yMat, k)
     [row , ~] = size(xMat);
     weights = zeros(row, row);
     for i = 1:1:row
         diffMat = point - xMat(i, :);
         weights(i,i) = exp(diffMat * diffMat.' / (-2.0 * (k ^ 2)));    %计算权重对角矩阵
     end
     xTx = xMat.' * (weights * xMat);   % 奇异矩阵不能计算
     if det(xTx) == 0
         printf('This matrix is singular, cannot do inverse');
     end
     theta = xTx^-1 * (xMat.' * (weights * yMat));          %计算回归系数
     yHatPoint = point * theta;
end

(2) 导入数据并绘图

clc;
clear all;
data = load('D:\\ex0.txt');
x = data(:,1:2);  
y = data(:,3);
y_hat = lwlr(x, x, y, 0.01);    % 调用局部加权线性回归函数计算估计值
x_axis = (x(:, 2))';            % x坐标轴为x矩阵的第二列数据
y_axis = y';                
sort(x_axis);                   % 对x向量升序排列
plot(x_axis, y, 'r.');          % 散点图
hold on
plot(x_axis, y_hat, 'b.');      %回归曲线图

(3) 分别选取k=1.0，0.01，0.003绘图如下：

图1 k=1.0 欠拟合

图2 k=0.01 最佳拟合

图3 k=0.003 过拟合

2.Python版
(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 i in range(m):
        diffMat = testPoint - xMat[i, :]
        weights[i, i] = 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
    theta = xTx.I * (xMat.T * (weights * yMat))                     # 计算回归系数
    return testPoint * theta

# 对所有点计算估计值
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

(2)导入数据并绘图

import matplotlib.pyplot as plt
from locally_weighted_linear_regression import LocallyWeightedLinearRegression as lwlr
from numpy import *

xArr, yArr = lwlr.loadDataSet('ex0.txt')
yHat = lwlr.lwlrTest(xArr, xArr, yArr, 0.01)
xMat = mat(xArr)
strInd = xMat[:, 1].argsort(0)
xSort = xMat[strInd][:, 0, :]

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

(3)分别选取k=1.0，0.01，0.003绘图如下：

图1 k=1.0 欠拟合

图2 k=0.01 最佳拟合

图3 k=0.003 过拟合

四、结果分析：很明显，当k=1.0时，回归曲线与散点图欠拟合(underfitting)，此时权重很大，如同将所有数据视为等权重，相当于普通线性回归模型；当k=0.01时得到了很好的效果；当k=0.003时，回归曲线与散点图过拟合(overfitting)，由于纳入了太多的噪声点，拟合的直线与数据点过于贴近，此时也不属于理想的模型。