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