本文大部分理论参考该博客【算法】局部加权回归（Lowess）仅用于自己学习记录
算法应用场景
Lowess局部加权回归算法本质作用就是拟合数据的趋势线，常见用于解决预测问题和平滑问题。
- 在解决预测问题时，利用趋势线来做预测数据，适用于周期性和波动性的数据；
- 在做数据平滑的时候，对有趋势或季节性的数据不能简单使用均值正负三倍标准差以外做异常值剔除，因此将趋势线作为基线，剔除偏离基线距离较远的真正异常值。
算法原理
核心思想
局部加权回归（Lowess）的大致思路是：以一个点x为中心，向前后截取一段长度为frac的数据，对于该数据用权值函数w做一个加权线性回归，记（x，y）为该回归线的中心值，其中y为拟合后曲线对应值。对于所有的n个数据点则可以做出n条加权回归线，每条回归线的中心值的连线则为这段数据的Lowess曲线。
权值函数选择
选择权值函数大致思路是：希望W(x)大于0，且作用域为[-1,1]，且为对称函数，该函数对于中间(0处)的值较大，两边(-1\1)处值较小。中间的权值较高，对于加权回归的影响较大；[-1,1]的原因是，对于任意不规则的数据段，可以压缩映射到[-1,1]，方便处理。
W
(
x
)
=
{
(
1
−
∣
x
∣
3
)
3
,
for
∣
x
∣
<
1
0
,
for
∣
x
∣
≥
1
W(x)= \begin{cases} (1-|x|^3)^3, &\text{for} |x|<1\\ 0, &\text{for} |x|≥1 \end{cases}
W(x)={(1−∣x∣3)3,0,for∣x∣<1for∣x∣≥1
B
(
x
)
=
{
(
1
−
∣
x
∣
2
)
2
,
for
∣
x
∣
<
1
0
,
for
∣
x
∣
≥
1
B(x) = \begin{cases} (1-|x|^2)^2, &\text{for} |x|<1\\ 0, &\text{for} |x|≥1 \end{cases}
B(x)={(1−∣x∣2)2,0,for∣x∣<1for∣x∣≥1
二次与三次函数的区别在于，三次函数对于周围权值降速更快，在平滑最初时候效果好，且适用于大多数分布，但增加了残差的方差。
对于权值函数选取，第一次迭代适用W函数（三次函数），之后迭代使用B函数（二次函数）。
权值函数使用
- 使用权值函数W(x)；
- 数据段[ d 1 , d 2 ] ，映射成[−1,1]对应的坐标；
- 带入函数W(x)，计算出每个点对应的w i ；
- 使用加权回归得出模型： Y ^ = X ( X T w X ) − 1 X T w Y \hat Y = X ( X ^T w X ) ^{-1}X ^T w Y Y^=X(XTwX)−1XTwY
残差计算与迭代
首先，原值为y ，预测值为y ^ ，残差为e = y − y ^，记s 为∣ e i ∣的中位数。鲁棒性的权值调整附加值δ k = W ( e k/ 6 s ) ，修正后的权值为δ k w k 。
迭代过程为：
1.使用W函数（三次函数）作为权值函数，求出w i 。
2.将w i 带入加权回归计算出y ^。
3.求出e = y − y 和s 。
4.以B函数作为修正权值函数，求出δ k = B ( e k/ 6 s） ，计算出δ k w k。
5.将δ k w k 作为修正权值，重复2、3、4步骤。
参数说明
- x_data,数据类型为列表，存储x的值
- y_data,数据类型为列表，存储y的值
- 长度frac，应该截取多长的作为局部处理，frac为原数据量的比例，为浮点数；
- 迭代次数iter，为int类型，函数基本会在3次左右收敛，一般取3
- 回归间隔delta，浮点数，为0时每个点都算一次加权回归；推荐当数据点N>5000时，delta=0.01*N，每隔delta个点进行一次加权回归，中间点采用：线性插值、二次插值、三次插值等方法。
调库简单使用
import statsmodels.api as sm
import matplotlib.pyplot as plt
import pandas as pd
data = pd.read_csv("Data.csv")
x_data = list(data.Date)
y_data = list(data.Value)
#直接调库使用lowess，输入都是一维数组，返回为点坐标的列表
lowess = sm.nonparametric.lowess(y_data, x_data, frac=0.1, it=3, delta=0.0)
x1 = [point[0] for point in lowess]
y1 = [point[1] for point in lowess]
plt.plot(x1,y1)
调库的缺点
- 权值函数不可调
- 使用delta时插值函数不可调节
前置准备：
- 编写截取数据函数get_points（输入x列表，y列表，中心x，截取比例frac；返回截取后的x列表，y列表）
- 确定权值函数W(x)和B(x)
- 编写加权回归函数weighted_linear_regression（x，y列表和权值列表w_list；返回 y ^ \hat{y} y^的列表）
- 编写计算新权值的函数new_W（输入 y ^ \hat{y} y^、y和权值列表w_list；输出残差、残差中位数和新权值列表）用于迭代过程
编程思路：
- 先检查输入的点坐标数据是否为顺序排列
- 循环遍历x
- 截取frac比例的数据作为列表
- 对列表内的点求权值函数并且做加权回归返回出 y ^ \hat{y} y^
- 进入迭代循环计算残差e、残差的中位数s，利用权值函数B(x)计算权值修正附加值，并输出新的权值
- 再计算得出新的 y ^ \hat{y} y^，重复至迭代结束（一般只迭代3次）
- 返回结果的中心点
python实现
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
####截取数据函数:以x为中心，向前后截取一段长度为frac的数据
def get_points(x_data,y_data,x,frac):
num = len(x_data)*frac
index = x_data.index(x)
if int(index-num//2) >= 0 and int(index+num//2) <= len(x_data):
x_list = x_data[int(index-num//2):int(index+num//2)]
y_list = y_data[int(index-num//2):int(index+num//2)]
elif int(index-num//2) < 0:
x_list = x_data[0:int(index+num//2)]
y_list = y_data[0:int(index+num//2)]
else:
x_list = x_data[int(index-num//2):len(x_data)]
y_list = y_data[int(index-num//2):len(x_data)]
return x_list,y_list
def weighted_linear_regression(x_list, y_list, w_list):
# 将输入列表转换为NumPy数组以进行数学操作
x = np.array(x_list)
y = np.array(y_list)
w = np.array(w_list)
# 计算加权线性回归的权重参数
# 这里使用了带权重的最小二乘法
#x要从列表转为设计矩阵（（（（（（（（（（（（（（（（（这里困扰了很久
X = np.vstack((x, np.ones_like(x))).T
#权重列表也需要转换为对角矩阵
W = np.diag(w)
'''
A = np.dot(X.T, np.dot(W, X))
B = np.dot(X.T, np.dot(W, y))
params = np.linalg.solve(A, B)
# 计算预测值
y_hat = np.dot(X, params)
'''
#也可直接写成一条
y_hat = X @ np.linalg.inv(X.T @ W @ X)@ X.T @ W @ y
return y_hat
####求残差、中位数、以及新权值
def new_W(y_hat,y_list,w_list):
e = y_list - y_hat ##计算对应的残差列表
e_sort = sorted(abs(e))
n = len(e)
##s为误差的中位数
if n % 2 == 1:
s = e_sort[(n - 1)//2]
else:
s = (e_sort[n//2]+e_sort[n//2+1])/2
new_w_list = np.array(B(e/(6*s)))*np.array(w_list)
return e,s,new_w_list
####二次权值函数B（x）
def B(x_list):
#用B函数求权值调整附加值时不需要映射到[-1,1]上，否则会使曲线不平滑
#x_list = [(x - min(x_list)) / (max(x_list) - min(x_list)) * 2 - 1 for x in x_list]
w_list = []
for i in x_list:
if abs(i) < 1:
w_list.append((1-i**2)**2)
else:
w_list.append(0)
return w_list
####三次权值函数W（x）
def W(x_list):
##先将x_list映射到[-1,1]上:对任意不规则数据段压缩映射到[-1,1]上不会改变原有序列的顺序
x_list = [(x - min(x_list)) / (max(x_list) - min(x_list)) * 2 - 1 for x in x_list]
w_list = []
for i in x_list:
if abs(i) < 1:
w_list.append((1 - abs(i) ** 3) ** 3)
else:
w_list.append(0)
return w_list
def mylowess(x_data,y_data,frac,iter):
result = []
#在x为中心截取代码段之前应该先检验x_data是否有序（后面都默认在x_data有序）
sorted_data = sorted(zip(x_data,y_data),key= lambda p: p[0])
x_data,y_data = zip(*sorted_data)
for x in x_data:
x_list, y_list = get_points(x_data, y_data, x, frac)
w_list = W(x_list)
y_hat = weighted_linear_regression(x_list, y_list, w_list)
for it in range(iter):
e, s, w_list = new_W(y_list, y_hat, w_list)
y_hat = weighted_linear_regression(x_list, y_list, w_list)
result.append(y_hat[x_list.index(x)])
return x_data, result
data = pd.read_csv("PolarArea_Data.csv")
x_data = list(data.Date)
y_data = list(data.Value)
frac = 0.1
iter = 3
plt.scatter(x_data,y_data)
for i in [1,2,3,4,5]:
x,y_fit = mylowess(x_data,y_data,frac,i)
plt.plot(x,y_fit, label = f"iter={i}")
plt.legend()
plt.show()
R语言实现
# 截取数据函数
get_points <- function(x_data, y_data, x, frac) {
num <- length(x_data) * frac
index <- which(x_data == x)
if (index - num / 2 >= 1 && index + num / 2 <= length(x_data)) {
x_list <- x_data[(index - num / 2):(index + num / 2)]
y_list <- y_data[(index - num / 2):(index + num / 2)]
} else if (index - num / 2 < 1) {
x_list <- x_data[1:(index + num / 2)]
y_list <- y_data[1:(index + num / 2)]
} else {
x_list <- x_data[(index - num / 2):length(x_data)]
y_list <- y_data[(index - num / 2):length(x_data)]
}
return(list(x_list, y_list))
}
# 加权线性回归函数
weighted_linear_regression <- function(x_list, y_list, w_list) {
X <- cbind(x_list, 1)
W <- diag(w_list)
params <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% y_list
y_hat <- X %*% params
return(y_hat)
}
# 二次权值函数B(x)
B <- function(x_list) {
w_list <- sapply(x_list, function(x) {
if (abs(x) < 1) {
return((1 - x^2)^2)
} else {
return(0)
}
})
return(w_list)
}
# 三次权值函数W(x)
W <- function(x_list) {
x_list <- (x_list - min(x_list)) / (max(x_list) - min(x_list)) * 2 - 1
w_list <- sapply(x_list, function(x) {
if (abs(x) < 1) {
return((1 - abs(x)^3)^3)
} else {
return(0)
}
})
return(w_list)
}
# 求残差、中位数、以及新权值
new_W <- function(y_hat, y_list, w_list) {
e <- y_list - y_hat
e_sort <- sort(abs(e))
n <- length(e)
if (n %% 2 == 1) {
s <- e_sort[(n - 1) %/% 2 + 1]
} else {
s <- (e_sort[n %/% 2] + e_sort[n %/% 2 + 1]) / 2
}
new_w_list <- B(e / (6 * s)) * w_list
return(list(e, s, new_w_list))
}
# mylowess函数
mylowess <- function(x_data, y_data, frac, iter) {
result <- numeric(length(x_data))
# 检查x_data是否有序
sorted_indices <- order(x_data)
x_data <- x_data[sorted_indices]
y_data <- y_data[sorted_indices]
for (x in x_data) {
points <- get_points(x_data, y_data, x, frac)
x_list <- points[[1]]
y_list <- points[[2]]
w_list <- W(x_list)
y_hat <- weighted_linear_regression(x_list, y_list, w_list)
for (it in 1:iter) {
new_w <- new_W(y_hat, y_list, w_list)
e <- new_w[[1]]
s <- new_w[[2]]
new_w_list <- new_w[[3]]
y_hat <- weighted_linear_regression(x_list, y_list, new_w_list)
}
result[which(x_data == x)] <- y_hat[which(x_list == x)]
}
return(result)
}
# 示例调用
data <- read.csv("PolarArea_Data.csv")
x_data <- data$Date
y_data <- data$Value
frac <- 0.1
iter <- 3
plot(x_data, y_data, pch = 19, col = "blue", xlab = "Date", ylab = "Value")
result <- mylowess(x_data, y_data, frac, iter)
lines(x_data, result, type = "l", col = "red")
# 添加标题、图例等其他绘图参数
title(main = "Scatter Plot with LOESS Fit", sub = "", xlab = "Date", ylab = "Value")
legend("topright", legend = c("Scatter", "LOESS Fit"), col = c("blue", "red"), pch = c(19, NA), lty = c(NA, 1), cex = 0.8)
结果比较
python调库使用迭代5次
mylowess迭代5次
可以看出mylowess还存在不小问题，每次迭代改变的幅度过小以至于几乎看不出有什么改变，还有很大的完善空间;另外，原调库函数中包含的间隔插值在这里没有实现，还有待补充。
反思与改进
关于列表映射到[-1,1]的步骤代码有误，应修改如下：
####三次权值函数W（x）
def W(xc,x_list):
##先将x_list映射到[-1,1]上:对任意不规则数据段压缩映射到[-1,1]上不会改变原有序列的顺序
x_list = [(x - xc) / (max(x_list) - min(x_list)) for x in x_list]
w_list = []
for i in x_list:
if abs(i) < 1:
w_list.append((1 - abs(i) ** 3) ** 3)
else:
w_list.append(0)
return w_list
运行结果
图像貌似更接近调库函数，但却不够平滑。由于数学水平有限，只知道当时错误乘2其实使得靠近中心点权重更大，远离中心点的权重更小
####三次权值函数W（x）
def W(xc,x_list):
##先将x_list映射到[-1,1]上:对任意不规则数据段压缩映射到[-1,1]上不会改变原有序列的顺序
x_list = [(x - xc) / (max(x_list) - min(x_list))*2 for x in x_list]
w_list = []
for i in x_list:
if abs(i) < 1:
w_list.append((1 - abs(i) ** 3) ** 3)
else:
w_list.append(0)
return w_list
还有一个问题就是对于稀疏点拟合曲线时，x值并不是连续的，按照比例取一定数量的点会出现问题，或许需要加一个距离筛选的判断
210
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
