【算法】局部加权回归(Lowess)

注:以下公式为latex公式,可看:http://blog.csdn.net/longgb123/article/details/79520982

一、简介

1.1 预测问题

对于预测问题,回归中最简单的线性回归,是以线性的方法拟合出数据的趋势。但是对于有周期性,波动性的数据,并不能简单以线性的方式拟合,否则模型会偏差较大,而局部加权回归(lowess)能较好的处理这种问题。可以拟合出一条符合整体趋势的线,进而做预测。

1.2 平滑问题

同时,局部加权回归(lowess)也能较好的解决平滑问题。在做数据平滑的时候,会有遇到有趋势或者季节性的数据,对于这样的数据,我们不能使用简单的均值正负3倍标准差以外做异常值剔除,需要考虑到趋势性等条件。使用局部加权回归,可以拟合一条趋势线,将该线作为基线,偏离基线距离较远的则是真正的异常值点。
实际上,局部加权回归(Lowess)主要还是处理平滑问题的多,因为预测问题,可以有更多模型做的更精确。但就平滑来说,Lowess很直观而且很有说服力。

二、算法讲解

2.1 算法思想

局部加权回归(Lowess)的大致思路是:以一个点$x$为中心,向前后截取一段长度为$frac$的数据,对于该段数据用权值函数$w$做一个加权的线性回归,记$(x,\hat{y})$为该回归线的中心值,其中$\hat{y}$为拟合后曲线对应值。对于所有的$n$个数据点则可以做出$n$条加权回归线,每条回归线的中心值的连线则为这段数据的Lowess曲线。

2.2 参数讲解

在这个思路中,能提取出的可调参数则是:
1.长度$frac$,应该截取多长的作为局部处理,$frac$为原数据量的比例;
2.权值函数$w$,使用什么样的权值函数$w$合适;
3.迭代次数$it$,在进行一次局部回归后,是否需要迭代,再次做回归;
4.$delta$回归间隔,是否真的每个点都需要算一次加权回归,能否隔$delta$距离算一次,中间没算的用插值替换即可。

在了解了算法算法的大致思想和可调参数以后,你可以马上上手使用statsmodels.api.nonparametric.lowess了。使用方法如下:

import statsmodels.api as sm
lowess = sm.nonparametric.lowess
result = lowess(y, x, frac=0.2, it=3, delta=0.0)

但是,在statsmodels中,你会发现:1、权值w函数你是不可调的;2、在用了$delta$之后,插值函数你是不可调的。
总之就是,黑盒子,很多你都不可调的。而且它还有一个非常严重的问题,具体可看本文3.3效果对比。
接下来就是看相关文档,了解思路,之后,你可以自己写一个lowess,而且速度不会慢。

相关文档,statsmodels就给出了比较权威的参考:《Cleveland, W.S. (1979) “Robust Locally Weighted Regression and Smoothing Scatterplots”. Journal of the American Statistical Association 74 (368): 829-836.》。
文章是《鲁棒性的加权回归》,即原始加权基础上迭代,增加鲁棒性。网上还有一些其他的lowess讲解,我看了,和这个不太一样,可以选择性阅读。

2.3 权值函数

理解了lowess之后,可以明白,其实权值函数并不是固定的,只要满足一定的规则条件即可(当然并也非强制),条件如下:
1.$W(x)>0 ; for ; |x|<1$;
2.$W(-x)=W(x)$;
3.$W(x) ; is ; a ; nonincreasing ; function ; for ; x \geqslant 0$;
4.$W(x)=0 ; for ; |x| \geqslant 1$

选择该类函数大致思路是:希望$W(x)$大于0,且作用域为[-1,1],且为对称函数,该函数对于中间(0处)的值较大,两边(-1\1)处值较小。
选择思路是,中间的权值较高,对于加权回归的影响较大;[-1,1]的原因是,对于任意不规则的数据段,可以压缩映射到[-1,1],方便处理。

权值函数如,B函数(二次函数):
$$
B(x)=\left{
\begin{aligned}
& (1-x{2}){2}, & for ; |x| < 1 \
& 0, & for ; |x| \geqslant 1
\end{aligned}
\right.
$$

W函数(三次函数):
$$
W(x)=\left{
\begin{aligned}
& (1-x{3}){3}, & for ; |x| < 1 \
& 0, & for ; |x| \geqslant 1
\end{aligned}
\right.
$$

二次与三次函数的区别在于,三次函数对于周围权值降速更快,在平滑最初时候效果好,且适用于大多数分布,但增加了残差的方差。
对于权值函数选取,第一次迭代适用W函数(三次函数),之后迭代使用B函数(二次函数)。
权值函数的使用:
1、使用权值函数$W(x)$;
2、数据段$[d_{1},d_{2}]$,映射成$[-1,1]$对应的坐标;
3、带入函数$W(x)$,计算出每个点对应的$w_{i}$;
4、使用加权回归得出模型:$\hat{Y}=X(X^{T}w X){-1}X{T}w Y$(推导见我的另一篇博客:线性回归,加权回归,推导过程

2.4 回归迭代

上面讲了权值函数的选取和使用,提到了迭代,这里讲解怎么迭代。
首先,原值为$y$,预测值为$\hat{y}$,残差为$e=y-\hat{y}$,记$s$为$|e_{i}|$的中位数。鲁棒性的权值调整附加值$\delta_{k}=W(\frac{e_{k}}{6s})$,修正后的权值为$\delta_{k}w_{k}$。
迭代过程为:
1.使用W函数(三次函数)作为权值函数,求出$w_{i}$。
2.将$w_{i}$带入加权回归计算出$\hat{y}$
3.求出$e=y-\hat{y}$和$s$
4.以B函数作为修正权值函数,求出$\delta_{k}=B(\frac{e_{k}}{6s})$,计算出$\delta_{k}w_{k}$
5.将$\delta_{k}w_{k}$作为修正权值,重复2、3、4步骤

该迭代没有明确的终止条件,据大量实验得知,原文中提到是2次迭代就基本收敛了,我做实验的时候,3次左右基本收敛,根文中描述差不多。

2.5 间隔回归,中间插值

在使用局部加权回归的时候,如果每个点都使用一次加权回归,则会比较耗时,所以有了,对于部分点使用加权回归,而未使用加权回归的点采用插值法处理,速度会增快很多,同时不会影响太大效果。
可以每间隔$delta$个点使用一次加权回归,中间点采用:线性插值、二次插值、三次插值等方法。
statsmodels推荐当数据点$N > 5000$的时候,选择$delta = 0.01 * N$
而我一般是设置$delta=4$,同时,我的插值函数是线性插值,因为我不希望插值出负数,大家可以根据自己需求选择。

2.6 其他参数

长度$frac$比例,文章给出的适用比例是:0.49,statsmodels默认的是:0.666。
同时,文章中还有一个参数是,使用的是多项式加权回归进行拟合,最高次是$d$,而进行讨论给出的结论是$d=1$
时候,即为线性回归,适合基本大多数场景。所以本博客一开始就使用线性加权回归介绍,感兴趣的可以去看看原文。

三、实验效果

3.1 效果

上面讲了整个思路,和详细的参数意义等,看了之后写出代码应该不难。这里作个伪代码总结:

data_x # x轴数据
data_y # y轴数据
map(x in data_x):               # 对每个x
    x_list = getRange(x, data_x, frac)  # 以x为中心,按照frac的比例截取数据
    w = calFuncW(x_list)                # 以w函数计算权值函数
    y_hat = weightRegression(x_list, data_y[x_list], w)         # 计算y_hat
    for it in iters:                    # 迭代iters次
        err, s, new_w = calNewWeight(y_hat, data_y[x_list], B)  # 使用B函数计算,迭代的new_w
        y_hat = weightRegression(x_list, data_y[x_list], new_w) # 使用new_w计算y_hat
    result = y_hat[x]                   # 取结果中心点

3.2 效率

上面自写的lowess,在spark上运行,executors15g4cores,计算30W8数据,每条数据90的数据条长度,10分钟即可完成(包括数据读写),个人觉得效率已经挺高的。

3.3 效果对比

statsmodels进行对比,右图为使用statsmodels的函数的结果,左图为自己根据文章写的运行结果,都是同一数据。

image.png

可以看出,对于一些特殊数据,直接使用statsmodels会有非常意外的结果,而且随着迭代进行,依旧不会收敛。

image.png

有的甚至会出现平滑结果很意外,如右图绿线最低部。

image.png

随着迭代进行,产生了更加剧烈的震荡。
所以,结论是:慎用statsmodelslowess函数。自己写的更可靠,而且自己可调部分更多。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,029评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,238评论 3 388
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 159,576评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,214评论 1 287
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,324评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,392评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,416评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,196评论 0 269
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,631评论 1 306
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,919评论 2 328
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,090评论 1 342
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,767评论 4 337
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,410评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,090评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,328评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,952评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,979评论 2 351

推荐阅读更多精彩内容