局部加权线性回归

最新推荐文章于 2024-02-14 00:00:00 发布

王亦傲

最新推荐文章于 2024-02-14 00:00:00 发布

阅读量1.9k

点赞数

本文链接：https://blog.csdn.net/weixin_44863781/article/details/100125575

版权

一、原理
局部加权线性回归是在一般线性回归的基础上进行改进，以解决一般线性回归的欠拟合问题。例如，如图所示
在这里插入图片描述
蓝点代表真实的样本，红线代表我们通过线性回归得到的模型，显然这个模型欠拟合了。那么，我们如何对对这个模型进行改进，已解决欠拟合问题呢？
我们知道对于一个样本点，与具有相同数据结构的数据在回归预测中预测值也相近。于是，我们可以让训练集中远离测试样本的训练样本对测试样本的影响减少，让训练集中接近测试样本的训练样本对测试样本的影响增加。那如何达成这个目标呢？我们知道一般线性回归的平方代价函数为 $J(θ)=\frac{1}{2m}\sum_{i=1}^m(h_θ(x_i)-y_i)^2$ ，我们给求和项中的每一项加上一个权重，即 $J(θ,x)=\frac{1}{2m}\sum_{i=1}^mW(x,x_i)(h_θ(x_i)-y_i)^2$ ,其中 $x$ 为待测样本。这个权重 $W(x,x_i)$ 具有这样的性质，待测样本 $x$ 与训练样本 $x_i$ 相近，则 $W(x,x_i)$ 大；否则， $W(x,x_i)$ 小。这里 $W(x,x_i)$ 的具体形式为：
$W(x,x_i)=exp(-\frac{|x_i-x|}{2k^2})$
其中参数 $k$ 调控着权重函数 $W(x,x_i)$ 的陡峭程度。如图所示为不同k值时权重函数 $W$ ,其中自变量为 $x_i-x|$ ,即 $x-x_i$ 的一维范数。
在这里插入图片描述
$k$ 越小，这个图形越陡峭，也就是说，当待测样本 $x$ 和训练样本 $x_i$ 相近会被分配很大的权重，而 $x$ 和 $x_i$ 相距较远会被分配的权重很小，几乎为0。 $k$ 较大时，图像相对较为均匀，则权重分配较为均匀。更显而易见的说法是， $k$ 控制着参与训练的样本的多少，当 $k$ 很小时，远离待测样本 $x$ 的训练样本被分配的权重几乎为0，也就是不参与训练，只有 $x$ 周围几个相近的样本参与训练。所以 $k$ 的选取尤为重要，如果 $k$ 值过小，只有少数样本参与训练必然造成过拟合;反之,这会导致欠拟合。
我们知道局部加权线性回归的代价函数为 $J(θ,x)=\frac{1}{2m}\sum_{i=1}^mW(x,x_i)(h_θ(x_i)-y_i)^2$ ，如果写成矩阵形式则有：
$J(θ)=(Xθ-y)^TW(Xθ-y)$
权重矩阵 $W$ 为一个m×m的矩阵，这是一个对角矩阵。其第i行第i列为样本 $x_i$ 的权重。
对 $J (θ)$ 关于 $θ$ 求导，得：
$\frac{∂J(θ)}{∂θ}=X^TW(Xθ-y)=0$
得到 $θ$ 的矩阵形式：
$θ=(X^TWX)^{-1}X^TWy$
最后，局部加权线性回归没有显性的训练过程，只有当给出测试样本时才能开始训练，属于懒惰学习。这一点和kNN算法一样。所以如果我们有一大堆新样本，会造成大量的训练过程。

二、实际应用
下面我们编写程序并用于一组数据。
①导入模块

import numpy as np
import matplotlib.pyplot as plt

②处理数据
首先查看原始数据：

file=open('F:/MachineLearning/data/lwlr.txt')
print(file.read())

得到：

1.000000        0.067732        3.176513
1.000000        0.427810        3.816464
1.000000        0.995731        4.550095
1.000000        0.738336        4.256571
1.000000        0.981083        4.560815
1.000000        0.526171        3.929515
1.000000        0.378887        3.526170
1.000000        0.033859        3.156393
1.000000        0.132791        3.110301
1.000000        0.138306        3.149813
1.000000        0.247809        3.476346
...

这里的第一列全为1，则对我们的回归不产生影响，我们可以去掉。但是由于常数项的存在我们还要在数据中加一列全是1的数组，所以第一列可以不用去掉。第三列实际是回归值，我们将它单独取出。
下面为数据处理程序：

def prepared_data(path):
    file=open(path)
    data_list=[]
    fr = file.readlines()
    for data in fr:
        data = data.split()
        data = [float(i) for i in data]
        data_list.append(data)
    data_list=np.array(data_list)
    y=data_list[:,-1]
    data_list=data_list[:,:-1]
    return data_list,y

第2—4行，读取数据并按行分隔。第5—8行，依次对每一行数据按间隔分隔，这里需要注意，我们列表里的数字是字符串类型，所以还要转化成数值类型。最后取出最后一行，并将数据类型从列表转化为numpy数组。
查看数据处理情况：

X,y=prepared_data('F:/MachineLearning/data/lwlr.txt')
print(X)
print(y)

得到：

[[1.       0.067732]
 [1.       0.42781 ]
 [1.       0.995731]
 [1.       0.738336]
 [1.       0.981083]
 [1.       0.526171]
 [1.       0.378887]
 [1.       0.033859]
 [1.       0.132791]
 [1.       0.138306]
 [1.       0.247809]
 [1.       0.64827 ]
 [1.       0.731209]
 [1.       0.236833]
 [1.       0.969788]
 [1.       0.607492]
 [1.       0.358622]
 [1.       0.147846]
 [1.       0.63782 ]
 [1.       0.230372]
 [1.       0.070237]
 [1.       0.067154]
 [1.       0.925577]
 [1.       0.717733]
 [1.       0.015371]
 [1.       0.33507 ]
 [1.       0.040486]
 [1.       0.212575]
 [1.       0.617218]
 [1.       0.541196]
 [1.       0.045353]
 ...
[3.176513 3.816464 4.550095 4.256571 4.560815 3.929515 3.52617  3.156393
 3.110301 3.149813 3.476346 4.119688 4.282233 3.486582 4.655492 3.965162
 3.5149   3.125947 4.094115 3.476039 3.21061  3.190612 4.631504 4.29589
 3.085028 3.44808  3.16744  3.364266 3.993482 3.891471 3.143259 3.114204
 3.851484 4.621899 4.580768 3.620992 3.580501 4.618706 3.676867 4.641845
 3.175939 4.26498  3.558448 3.436632 3.831052 3.182853 3.498906 3.946833
 3.900583 4.238522 4.23308  3.521557 3.203344 4.278105 3.555705 3.502661
 3.859776 4.275956 3.916191 3.587961 3.183004 4.225236 4.231083 4.240544
 3.222372 4.021445 3.567479 3.56258  4.262059 3.208813 3.169825 4.193949
 3.491678 4.533306 3.550108 4.636427 3.557078 3.552874 3.494159 3.206828
 3.195266 4.221292 4.413372 4.184347 3.742878 3.201878 4.648964 3.510117
 3.274434 3.579622 3.489244 4.237386 3.913749 3.22899  4.286286 4.628614
 3.239536 4.457997 3.513384 3.729674 3.834274 3.811155 3.598316 4.692514
 4.604859 3.864912 3.184236 3.500796 3.743365 3.622905 4.310796 3.583357
 3.901852 3.233521 3.105266 3.865544 4.628625 4.231213 3.791149 3.968271
 4.25391  3.19471  3.996503 3.904358 3.503976 4.557545 3.699876 4.613614
 3.140401 4.206717 3.969524 4.476096 3.136528 4.279071 3.200603 3.299012
 3.209873 3.632942 3.248361 3.995783 3.563262 3.649712 3.951845 3.145031
 3.181577 4.637087 3.404964 3.873188 4.633648 3.154768 4.623637 3.078132
 3.913596 3.221817 3.938071 3.880822 4.176436 4.648161 3.332312 4.240614
 4.532224 4.557105 4.610072 4.636569 4.229813 3.50086  4.245514 4.605182
 3.45434  3.180775 3.38082  4.56502  3.279973 4.554241 4.63352  4.281037
 3.844426 3.891601 3.849728 3.492215 4.592374 4.632025 3.75675  3.133555
 3.567919 4.363382 3.560165 4.564305 4.215055 4.174999 4.58664  3.960008
 3.529963 4.213412 3.908685 3.585821 4.374394 3.213817 3.952681 3.129283]

下面通过图来查看数据情况：

fig1,ax1 = plt.subplots(figsize=(9,6))
ax1.scatter(X[:,-1],y,s=10)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
plt.show()

得到：
在这里插入图片描述

③训练器函数部分
首先编写一个得到参数 $θ$ 的函数：

def lwlr(xTest,X,y,k):
    m,d = X.shape
    W=np.zeros((m,m))
    for i in range(m):
        W[i,i]=np.exp(np.linalg.norm(xTest-X[i,:])/(-2*k**2))
    XwX = X.T @ W @ X
    if np.linalg.det(XwX)==0:
        print('XwX can not inverse')
    theta = np.linalg.inv(XwX) @ ((X.T @ W) @ y)
    return theta

这个函数需要输入我们待预测点xTest,训练数据X和y以及参数k。这个函数前面先计算权重矩阵 $W$ ，然后判断 $X^TWX$ 是否可逆。最后利用公式 $θ=(X^TWX)^{-1}X^TWy$ 得到我们的参数。
下面通过参数得到预测值：

def lwlrPredict(xTest,omega):
    predict_y = xTest @ omega
    return predict_y

最后给定测试集，得到测试集的预测值：

def lwlrTest(Xtest,X,y,k):
    predict_yList = []
    m_test,d_test = Xtest.shape
    for i in range(m_test):
        omega=lwlr(Xtest[i,:],X,y,k)
        predict_y = lwlrPredict(Xtest[i,:],omega)
        predict_yList.append(predict_y)
    return np.array(predict_yList)

④应用于数据
这里我们直接把训练集当做待测样本集，查看拟合情况。

predict_y=lwlrTest(X,X,y,0.08)
X_index=X[:,1].argsort()
X_sort=X[:,1][X_index]
predict_y=predict_y[X_index]
fig1,ax1 = plt.subplots(figsize=(9,6))
ax1.scatter(X[:,1],y,s=10)
ax1.plot(X_sort,predict_y,c='red')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
plt.show()

这里我们把参数 $k$ 设置为0.08,得到：
在这里插入图片描述
可以看到拟合明显较好。如果我们把参数变小设置为0.03，得到：

在这里插入图片描述
明显有些过拟合。我们再把参数调大些，调成1：

很明显又欠拟合了。所以对于参数 $k$ 的选取尤为重要。

三、补充：用梯度下降法实现
上面我们我们用了正则方程的方法实现了局部加权线性回归，其实还可以用梯度下降法的迭代方法去实现。
下面为这一方法的代码部分：

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as pyplot

def computeWeight(xTest,X,k):
    m,d=X.shape
    W=np.zeros((m,m))
    for i in range(m):
        W[i,i]=np.exp(-np.linalg.norm(xTest-X[i,:])/(2*k**2))
    return W

def cost(theta,xTest,X,y,k):
    m,d = X.shape
    W = computeWeight(xTest,X,k)
    cost = (X @ theta - y).T @ W @ (X @ theta - y) / (2 * m)
    return cost

def gradient(theta,xTest,X,y,k):
    m,d = X.shape
    W = computeWeight(xTest,X,k)
    gradient = X.T @ W @ (X @ theta - y) / m
    return gradient


def TestGradient(Xtest,X,y,k):
    mTest,dTest = Xtest.shape
    theta = np.zeros(dTest)
    predict_y = np.zeros(mTest)
    for j in range(mTest):
        res = opt.minimize(fun=cost, x0=theta, args=(Xtest[j,:],X, y,k), method='Newton-CG', jac=gradient)
        predict_y[j] = Xtest[j,:] @ res.x
    return predict_y

下面我们看看拟合情况：

from lwlr import prepared_data

X,y=prepared_data('F:/MachineLearning/data/lwlr.txt')


predict_y=TestGradient(X,X,y,0.08)

X_index=X[:,1].argsort()
X_sort=X[:,1][X_index]
predict_y=predict_y[X_index]
fig1,ax1 = plt.subplots(figsize=(9,6))
ax1.scatter(X[:,1],y,s=10)
ax1.plot(X_sort,predict_y,c='red')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
plt.show()

于是得到，参数 $k$ 设置为0.08时的拟合情况：
在这里插入图片描述
$k$ 为0.03时：

$k$ 为1时：

王亦傲

关注

0
点赞
踩
12

收藏

觉得还不错? 一键收藏
2
评论
局部加权线性回归

一、原理局部加权线性回归是在一般线性回归的基础上进行改进，以解决一般线性回归的欠拟合问题。例如，如图所示蓝点代表真实的样本，红线代表我们通过线性回归得到的模型，显然这个模型欠拟合了。那么，我们如何对对这个模型进行改进，已解决欠拟合问题呢？我们知道对于一个样本点，与具有相同数据结构的数据在回归预测中预测值也相近。于是，我们可以让训练集中远离测试样本的训练样本对测试样本的影响减少，让训练集中接...
复制链接

扫一扫